A systematic global investigation of differential charge radii has been performed within the CDFT framework for the first time. Theoretical results obtained with conventional covariant energy density functionals and separable pairing interaction are compared with experimental differential charge radii in the regions of the nuclear chart in which available experimental data crosses neutron shell closures at N = 28, 50, 82 and 126. The analysis of absolute differential radii of different isotopic chains and their relative properties indicate clearly that such properties are reasonably well described in model calculations in the cases when the mean-field approximation is justified. However, while the observed clusterization of differential charge radii of different isotopic chains is well described above the N=50 and N=126 shell closures, it is more difficult to reproduce it above the N=28 and N=82 shell closures because of possible deficiencies in underlying single-particle structure. The impact of the latter has been evaluated for spherical shapes and it was shown that the relative energies of the single-particle states and the patterns of their occupation with increasing neutron number have an appreciable impact on the evolution of the differential charge radii. It is shown that the kinks in the charge radii at neutron shell closures are due to the underlying single-particle structure and due to weakening or collapse of pairing at these closures. It is usually assumed that pairing is a dominant contributor to odd-even staggering (OES) in charge radii. Our analysis paints a more complicated picture. It suggests a new mechanism in which the fragmentation of the single-particle content of the ground state in odd-mass nuclei due to particle-vibration coupling provides a significant contribution to OES in charge radii.