  We present the first fully consistent calculation of inclusive $\pi^0$
production at forward rapidities at next-to-leading order (NLO) accuracy in
proton-lead collisions at $\sqrt{s}=8.16$ TeV within the Color Glass Condensate
approach. The center-of-mass energy dependence is determined by the
Balitsky-Kovchegov equation with the initial condition constrained at NLO
accuracy by the proton structure function measurements. We find that the LHCb
data further constrain this non-perturbative input, and that both the NLO
corrections to the evolution equation and to the impact factor have significant
effects on the nuclear modification factor $R_{p\mathrm{Pb}}$. The complete NLO
calculation is shown to be in qualitative agreement with the LHCb data, and the
predictions are found to be sensitive to the details of the applied deep
inelastic scattering (DIS) fit. These results demonstrate the need to perform
global analyses of DIS and LHC data in order to probe gluon saturation
phenomena at precision level at very small momentum fraction $x$.
