Splitting functions govern the scale evolution of parton distribution functions. Through a Mellin transformation, they are related to anomalous dimensions of twist-two operators in the operator product expansion. We study off-shell operator matrix element, where the physical operators mix under renormalization with other gauge-variant operators of the same quantum numbers. We devise a new method to systematically extract the Feynman rules resulting from those operators without knowing the operators themselves. As a first application of the new approach, we independently reproduce the well-known three-loop singlet splitting functions obtained from computations of on-shell quantities.