For an arbitrary renormalizable N = 1 supersymmetric theory with a single gauge coupling regularized by higher covariant derivatives we calculate the three-loop contribution to the anomalous dimension of the matter superfields proportional to the sixth powers of Yukawa couplings for a wide class of renormalization prescriptions. Next, taking into account that in the case of using the higher covariant derivatives supplemented by minimal subtractions of logarithms the NSVZ relation is valid in all orders, we calculate the corresponding four-loop contribution to the gauge β-function and investigate its scheme dependence. In particular, we construct the so-called “minimal” renormalization prescription for which the renormalization group functions have the simplest possible form provided that the exact equations for the gauge and Yukawa β-functions following from supersymmetry remain valid.