  The leading short-distance contributions to hadronic hard-scattering cross
sections in the operator product expansion are described by twist-two quark and
gluon operators. The anomalous dimensions of these operators determine the
splitting functions that govern the scale evolution of parton distribution
functions. In massless QCD, these anomalous dimensions can be determined
through the calculation of off-shell operator matrix elements, typically
performed in a covariant gauge, where the physical operators mix with
gauge-variant operators of the same quantum numbers. We derive a new method to
systematically extract the counterterm Feynman rules resulting from these
gauge-variant operators. As a first application of the new method, we rederive
the unpolarized three-loop singlet anomalous dimensions, independently
confirming previous results obtained with other methods. Employing a general
covariant gauge, we observe the explicit cancellation of the gauge parameter
dependence in these results.
