We present the MSHT20qed_an3lo parton distribution functions (PDFs). These result from the first global PDF analysis to combine QED and approximate N^3N3LO (aN^3N3LO) QCD corrections in the theoretical calculation of the PDF evolution and cross sections entering the fit. We examine the PDF impact, and find that the effect of QED is relatively mild in comparison to the aN^3N3LO corrections, although it should still be accounted for at the level of precision now required. These QED corrections are in addition found to roughly factorise from the QCD corrections; that is, their relative impact on the PDFs is roughly the same at NNLO and aN^3N3LO. The fit quality exhibits a very small deterioration at aN^3N3LO upon the inclusion of QED corrections, which is rather smaller than the deterioration observed at NNLO in QCD. The impact on several cross-sections at N^3N3LO is also examined, including the Higgs cross sections at N^3N3LO. Finally, a LO in QCD fit that includes QED corrections is also presented: the MSHT20qed_an3lo set.