  We calculate the unpolarized and polarized two--loop massless off--shell
operator matrix elements in QCD to $O(\varepsilon)$ in the dimensional
parameter in an automated way. Here we use the method of arbitrary high Mellin
moments and difference ring theory, based on integration-by-parts relations.
This method also constitutes one way to compute the QCD anomalous dimensions.
The presented higher order contributions to these operator matrix elements
occur as building blocks in the corresponding higher order calculations up to
four--loop order. All contributing quantities can be expressed in terms of
harmonic sums in Mellin--$N$ space or by harmonic polylogarithms in $z$--space.
We also perform comparisons to the literature.
