We present an updated global analysis of collinearly factorized nuclear parton distribution functions (PDFs) at next-to-leading order in perturbative QCD. In comparison to our previous fit, EPPS16, the present analysis includes more data from proton-lead collisions at the Large Hadron Collider: 5 TeV double-differential CMS dijet and LHCb D-meson data, as well as 8 TeV CMS W± data. These new data lead to significantly better-constrained gluon distributions at small and intermediate values of the momentum fraction x, confirming the presence of shadowing and antishadowing for gluons in large nuclei. In addition, we include Jefferson Lab measurements of deeply inelastic scattering which probe nuclear PDFs at large x and low virtualities. For the first time within the Hessian framework, we now also explore the uncertainties of nuclear PDFs due to the errors in the baseline proton PDFs. We release the results of our analysis as a new public parametrization of nuclear PDFs called EPPS21.