  We show that integrable structure of chiral random matrix models
incorporating global symmetries of QCD Dirac operators (labeled by the Dyson
index beta=1,2, and 4) leads to emergence of a connection relation between the
spectral statistics of massive and massless Dirac operators. This novel link
established for beta-fold degenerate massive fermions is used to explicitly
derive (and prove the random matrix universality of) statistics of low--lying
eigenvalues of QCD Dirac operators in the presence of SU(2) massive fermions in
the fundamental representation (beta=1) and SU(N_c >= 2) massive adjoint
fermions (beta=4). Comparison with available lattice data for SU(2) dynamical
staggered fermions reveals a good agreement.
