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 land prove the random matrix universality of) statistics of low-lying spectra of QCD Dirac operators in the presence of SU(2) massive fermions in the fundamental representation (beta = 1) and SU(N-c greater than or equal to 2) massive adjoint fermions (beta = 4). Comparison with available lattice data for SU(2) dynamical staggered fermions reveals a good agreement.