  We investigate the NSVZ relations for ${\cal 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 $\beta$-functions. For ${\cal 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 $\beta$-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
$\beta$-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.
