An atomic-orbital (AO) based formulation for calculating nuclear magnetic resonance chemical shieldings at the second-order Møller-Plesset perturbation theory level is introduced, which provides a basis for reducing the scaling of the computational effort with the molecular size from the fifth power to linear and for a specific nucleus to sublinear. The latter sublinear scaling in the rate-determining steps becomes possible by avoiding global perturbations with respect to the magnetic field and by solving for quantities that involve the local nuclear magnetic spin perturbation instead. For avoiding the calculation of the second-order perturbed density matrix, we extend our AO-based reformulation of the Z-vector method within a density matrix-based scheme. Our pilot implementation illustrates the fast convergence with respect to the required number of Laplace points and the asymptotic scaling behavior in the rate-determining steps.