Using Nf=2+1+1 lattice QCD, we determine the fermionic connected contributions to the first and second moments of the pion parton distribution function. Based on gauge configurations from the European Twisted Mass Collaboration, chiral and continuum extrapolations are performed using pion masses in the range of 230 to 500 MeV and three values of the lattice spacing. Finite volume effects are investigated using different volumes. In order to avoid mixing under renormalization for the second moment, we use an operator with two nonzero spatial components of momentum. Momenta are injected using twisted boundary conditions. Our final values read ⟨x⟩Rphys=0.2075(106) and ⟨x2⟩Rphys=0.163(33), determined at 2 GeV in the MS¯ scheme and with systematic and statistical uncertainties summed in quadrature.