We present results of our recent EPPS16 global analysis of NLO nuclear parton distribution functions (nPDFs). For the first time, dijet and heavy gauge boson production data from LHC proton– lead collisions have been included in a global fit. Especially, the CMS dijets play an important role in constraining the nuclear effects in gluon distributions. With the inclusion of also neutrino– nucleus deeply-inelastic scattering and pion–nucleus Drell–Yan data and a proper treatment of isospin-corrected data, we were able to free the flavor dependence of the valence and sea quark nuclear modifications for the first time. This gives us less biased, yet larger, flavor by flavor uncertainty estimates. The EPPS16 analysis indicates no tension between the data sets used, which supports the validity of collinear factorization and universal nPDFs for nuclear hard-collision processes in the kinematical range studied.