In pursuing the essential elements of nuclear binding, we compute ground-state properties of atomic nuclei with up to A=20 nucleons, using as input a leading order pionless effective field theory Hamiltonian. A variational Monte Carlo method based on a new, highly-expressive, neural-network quantum state ansatz is employed to solve the many-body Schrödinger equation in a systematically improvable fashion. In addition to binding energies and charge radii, we accurately evaluate the magnetic moments of these nuclei, as they reveal the self-emergence of the shell structure, which is not a priori encoded in the neural-network ansatz. To this aim, we introduce a novel computational protocol based on adding an external magnetic field to the nuclear Hamiltonian, which allows the neural network to learn the preferred polarization of the nucleus within the given magnetic field.