  We study the relevance of experimental data on heavy-flavor [$D^0$, $J/\psi$,
$B\rightarrow 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 $x\simeq 7\times 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.
