Analytic calculation of nuclear magnetic resonance chemical shielding tensors, based on gauge-including atomic orbitals, is implemented for double-hybrid density functional theory (DHDFT), using the resolution of the identity (RI) approximation for its second order Møller-Plesset perturbation theory (MP2) correlation contributions. A benchmark set of 15 small molecules, containing 1H, 13C, 15N, 17O, 19F, and 31P nuclei, is used to assess the accuracy of the results in comparison to coupled cluster and empirical equilibrium reference data, as well as to calculations with MP2, Hartree-Fock, and commonly used pure and hybrid density functionals. Investigated are also errors due to basis set incompleteness, the frozen core approximation, different auxiliary basis sets for the RI approximation, and grids used for the chain-of-spheres exchange integral evaluation. The DSD-PBEP86 double-hybrid functional shows the smallest deviations from the reference data with mean absolute relative error in chemical shifts of 1.9%. This is significantly better than MP2 (4.1%), spin-component-scaled MP2 (3.9%), or the best conventional density functional tested, M06L (5.4%). A protocol (basis sets, grid sizes, etc.) for the efficient and accurate calculation of chemical shifts at the DHDFT level is proposed and shown to be routinely applicable to systems of 100-400 electrons, requiring computation times 1-2 orders of magnitude longer than for equivalent calculations with conventional (pure or hybrid) density functionals.