We extend the existing leading (LO), next-to-leading (NLO), and next-to-next-to-leading order (NNLO) NNPDF4.0 sets of parton distribution functions (PDFs) to approximate next-to-next-to-next-to-leading order (aN3LO). We construct an approximation to the N3LO splitting functions that includes all available partial information from both fixed-order computations and from small and large x resummation, and estimate the uncertainty on this approximation by varying the set of basis functions used to construct the approximation. We include known N3LO corrections to deep-inelastic scattering structure functions and extend the FONLL general-mass scheme to Oαs3 accuracy. We determine a set of aN3LO PDFs by accounting both for the uncertainty on splitting functions due to the incomplete knowledge of N3LO terms, and to the uncertainty related to missing higher corrections (MHOU), estimated by scale variation, through a theory covariance matrix formalism. We assess the perturbative stability of the resulting PDFs, we study the impact of MHOUs on them, and we compare our results to the aN3LO PDFs from the MSHT group. We examine the phenomenological impact of aN3LO corrections on parton luminosities at the LHC, and give a first assessment of the impact of aN3LO PDFs on the Higgs and Drell–Yan total production cross-sections. We find that the aN3LO NNPDF4.0 PDFs are consistent within uncertainties with their NNLO counterparts, that they improve the description of the global dataset and the perturbative convergence of Higgs and Drell–Yan cross-sections, and that MHOUs on PDFs decrease substantially with the increase of perturbative order.