  We scrutinize the recent LHCb data for D$^0$-meson production in p+Pb
collisions within a next-to-leading order QCD framework. Our calculations are
performed in the SACOT-$m_{\rm T}$ variant of the general-mass
variable-flavour-number scheme (GM-VFNS), which has previously been shown to
provide a realistic description of the LHC p+p data. Using the EPPS16 and
nCTEQ15 nuclear parton distribution functions (PDFs) we show that a very good
agreement is obtained also in the p+Pb case both for cross sections and nuclear
modification ratios in the wide rapidity range covered by the LHCb data.
Encouraged by the good correspondence, we quantify the impact of these data on
the nuclear PDFs by the Hessian reweighting technique. We find compelling
direct evidence of gluon shadowing at small momentum fractions $x$, with no
signs of parton dynamics beyond the collinear factorization. We also compare
our theoretical framework to a fixed-order calculation supplemented with a
parton shower. While the two frameworks differ in the absolute cross sections,
these differences largely cancel in the nuclear modification ratios. Thus, the
constraints for nuclear PDFs appear solid.
