Precise quantification of the structure of nucleons is one of the crucial aims of hadronic physics for the coming years. The expected progress related to ongoing and planned experiments should be accompanied by calculations of partonic distributions from lattice QCD. While key insights from the lattice are expected to come for distributions that are difficult to access experimentally, it is important that lattice QCD can reproduce the well-known unpolarized parton distribution functions (PDFs) with full control over systematic uncertainties. One of the novel methods for accessing the partonic x-dependence is the pseudo-distribution approach, which employs matrix elements of a spatially-extended nonlocal Wilson-line operator of length z. In this paper, we address the issue of discretization effects, related to the necessarily nonzero value of the lattice spacing a, which start at first order in a as a result of the nonlocal operator. We use twisted mass fermions simulated at three values of the lattice spacing, at a pion mass of 370 MeV, and extract the continuum limit of isovector unpolarized PDFs. We also test, for the first time in the pseudo-distribution approach, the effects of the recently derived two-loop matching. Finally, we address the issue of the reliability of the extraction with respect to the maximal value of z.