A spherical Gaussian nuclear charge distribution model has been implemented for spin-free (scalar) and two-component (spin-orbit) relativistic density functional calculations of indirect NMR nuclear spin-spin coupling (J-coupling) constants. The finite nuclear volume effects on the hyperfine integrals are quite pronounced and as a consequence they noticeably alter coupling constants involving heavy NMR nuclei such as W, Pt, Hg, Tl, and Pb. Typically, the isotropic J-couplings are reduced in magnitude by about 10 to 15 % for couplings between one of the heaviest NMR nuclei and a light atomic ligand, and even more so for couplings between two heavy atoms. For a subset of the systems studied, viz. the Hg atom, Hg(2) (2+), and Tl--X where X=Br, I, the basis set convergence of the hyperfine integrals and the coupling constants was monitored. For the Hg atom, numerical and basis set calculations of the electron density and the 1s and 6s orbital hyperfine integrals are directly compared. The coupling anisotropies of TlBr and TlI increase by about 2 % due to finite-nucleus effects.