We investigate the NSVZ relations for $$ \mathcal{N} $$
N
= 1 supersymmetric gauge theories with multiple gauge couplings. As examples, we consider MSSM and the flipped SU(5) model, for which they easily reproduce the results for the two-loop β-functions. For $$ \mathcal{N} $$
N
= 1 SQCD interacting with the Abelian gauge superfield we demonstrate that the NSVZ-like equation for the Adler D-function follows from the NSVZ relations. Also we derive all-loop equations describing how the NSVZ equations for theories with multiple gauge couplings change under finite renormalizations. They allow describing a continuous set of NSVZ schemes in which the exact NSVZ β-functions are valid for all gauge coupling constants. Very likely, this class includes the HD+MSL scheme, which is obtained if a theory is regularized by Higher covariant Derivatives and divergences are removed by Minimal Subtractions of Logarithms. That is why we also discuss how one can construct the higher derivative regularization for theories with multiple gauge couplings. Presumably, this regularization allows to derive the NSVZ equations for such theories in all loops. In this paper we make the first step of this derivation, namely, the NSVZ equations for theories with multiple gauge couplings are rewritten in a new form which relates the β-functions to the anomalous dimensions of the quantum gauge superfields, of the Faddeev-Popov ghosts, and of the matter superfields. The equivalence of this new form to the original NSVZ relations follows from the extension of the non-renormalization theorem for the triple gauge-ghost vertices, which is also derived in this paper.