We present a determination of the parton distribution functions of the proton in which NLO and NNLO fixed-order calculations are supplemented by NLLx small-x resummation. Deep inelastic structure functions are computed consistently at NLO+NLLx or NNLO+NLLx, while for hadronic processes small-x resummation is included only in the PDF evolution, with kinematic cuts introduced to ensure the fitted data lie in a region where the fixed-order calculation of the hard cross-sections is reliable. In all other respects, the fits use the same methodology and are based on the same global dataset as the recent NNPDF3.1 analysis. We demonstrate that the inclusion of small-x resummation leads to a quantitative improvement in the perturbative description of the HERA inclusive and charm-production reduced cross-sections in the small x region. The impact of the resummation in our fits is greater at NNLO than at NLO, because fixed-order calculations have a perturbative instability at small x due to large logarithms that can be cured by resummation. We explore the phenomenological implications of PDF sets with small-x resummation for the longitudinal structure function $F_L$ at HERA, for parton luminosities and LHC benchmark cross-sections, for ultra-high energy neutrino-nucleus cross-sections, and for future high-energy lepton-proton colliders such as the LHeC.