We present a joint analysis of the four most prominent sodium-sensitive features (Na D, Na I λ8190Å, Na I λ1.14 μm, and Na I λ2.21 μm), in the optical and near-infrared spectral ranges, of two nearby, massive (σ ˜ 300 km s-1), early-type galaxies (named XSG1 and XSG2). Our analysis relies on deep Very Large Telescope/X-Shooter long-slit spectra, along with newly developed stellar population models, allowing for [Na/Fe] variations, up to ˜1.2 dex, over a wide range of age, total metallicity, and initial mass function (IMF) slope. The new models show that the response of the Na-dependent spectral indices to [Na/Fe] is stronger when the IMF is bottom heavier. For the first time, we are able to match all four Na features in the central regions of massive early-type galaxies finding an overabundance of [Na/Fe] in the range 0.5-0.7 dex and a bottom-heavy IMF. Therefore, individual abundance variations cannot be fully responsible for the trends of gravity-sensitive indices, strengthening the case towards a non-universal IMF. Given current limitations of theoretical atmosphere models, our [Na/Fe] estimates should be taken as upper limits. For XSG1, where line strengths are measured out to ˜0.8 Re, the radial trend of [Na/Fe] is similar to [α/Fe] and [C/Fe], being constant out to ˜0.5 Re, and decreasing by ˜0.2-0.3 dex at ˜0.8 Re, without any clear correlation with local metallicity. Such a result seems to be in contrast to the predicted increase of Na nucleosynthetic yields from asymptotic giant branch stars and Type II supernovae. For XSG1, the Na-inferred IMF radial profile is consistent, within the errors, with that derived from TiO features and the Wing-Ford band presented in a recent paper.