  Explicit diagrammatic calculation of evolution equations for high-twist
correlation functions is a challenge already at one-loop order in QCD coupling.
The main complication being quite involved mixing pattern of the so-called
non-quasipartonic operators. Recently, this task was completed in the
literature for twist-four nonsinglet sector. Presently, we elaborate on a
particular component of renormalization corresponding to the mixing of
gauge-invariant operators with QCD equations of motion. These provide an
intrinsic contribution to evolution equations yielding total result in
agreement with earlier computations that bypassed explicit analysis of Feynman
graphs.
