We study the relevance of experimental data on heavy-flavor [$D^0$, $J/psi$, $Brightarrow J/psi$ and $Upsilon(1S)$ mesons] production in proton-lead collisions at the LHC to improve our knowledge of the gluon-momentum distribution inside heavy nuclei. We observe that the nuclear effects encoded in both most recent global fits of nuclear parton densities at next-to-leading order (nCTEQ15 and EPPS16) provide a good overall description of the LHC data. We interpret this as a hint that these are the dominant ones. In turn, we perform a Bayesian-reweighting analysis for each particle data sample which shows that each of the existing heavy-quark(onium) data set clearly points --with a minimal statistical significance of 7 $sigma$-- to a shadowed gluon distribution at small $x$ in the lead. Moreover, our analysis corroborates the existence of gluon antishadowing. Overall, the inclusion of such heavy-flavor data in a global fit would significantly reduce the uncertainty on the gluon density down to $xsimeq 7times 10^{-6}$ --where no other data exist-- while keeping an agreement with the other data of the global fits. Our study accounts for the factorization-scale uncertainties which dominate for the charm(onium) sector.