Massive Early-Type Galaxies (ETGs) in the local Universe are believed to be the most mature stage of galaxy evolution. Their stellar population content reveals the evolutionary history of these galaxies. However, while state-of-the-art Stellar Population Synthesis (SPS) models provide an accurate description of observed galaxy spectra in the optical range, the modelling in the Near-Infrared (NIR) is still in its infancy. Here, we focus on NIR CO absorption features to show, in a systematic and comprehensive manner, that for massive ETGs, all CO indices, from H through to K band, are significantly stronger than currently predicted by SPS models. We explore and discuss several possible explanations of this 'CO mismatch', including the effect of intermediate-age, asymptotic-giant-branch-dominated, stellar populations, high-metallicity populations, non-solar abundance ratios, and the initial mass function. While none of these effects is able to reconcile models and observations, we show that ad hoc 'empirical' corrections, taking into account the effect of CO-strong giant stars in the low-temperature regime, provide model predictions that are closer to the observations. Our analysis points to the effect of carbon abundance as the most likely explanation of NIR CO line-strengths, indicating possible routes for improving the SPS models in the NIR.