We apply the Kurland-McGarvey (J. Magn. Reson. 1970, 2, 286) theory for the NMR shielding of paramagnetic molecules, particularly its special case limited to the ground-state multiplet characterized by zero-field splitting (ZFS) interaction of the form S·D·S. The correct formulation for this problem was recently presented by Soncini and Van den Heuvel (J. Chem. Phys. 2013, 138, 054113). With the effective electron spin quantum number S, the theory involves 2S+1 states, of which all but one are low-lying excited states, between which magnetic couplings take place by Zeeman and hyperfine interactions. We investigate these couplings as a function of temperature, focusing on both the high- and low-temperature behaviors. As has been seen in work by others, the full treatment of magnetic couplings is crucial for a realistic description of the temperature behavior of NMR shielding up to normal measurement temperatures. At high temperatures, depending on the magnitude of ZFS, the effect of magnetic couplings diminishes, and the Zeeman and hyperfine interactions become effectively averaged in the thermally occupied states of the multiplet. At still higher temperatures, the ZFS may be omitted altogether, and the shielding properties may be evaluated using a doublet-like formula, with all the 2S+1 states becoming effectively degenerate at the limit of vanishing magnetic field. We demonstrate these features using first-principles calculations of Ni(II), Co(II), Cr(II), and Cr(III) complexes, which have ZFS of different sizes and signs. A non-monotonic inverse temperature dependence of the hyperfine shift is predicted for axially symmetric integer-spin systems with a positive D parameter of ZFS. This is due to the magnetic coupling terms that are proportional to kT at low temperatures, canceling the Curie-type 1/kT prefactor of the hyperfine shielding in this case.