Solving the time-dependent Schrödinger equation (TDSE) for large molecular systems is a complicated task due to the inherent exponential scaling of the problem. One of the most successful and versatile methods for obtaining numerically converged solutions for small to medium-sized systems is multiconfiguration time-dependent Hartree (MCTDH). In a recent publication [J. Chem. Phys. 152, 084101 (2020)] we introduced a hierarchy of approximations to the MCTDH method which mitigate the exponential scaling by truncating the configuration space based on a maximum excitation level w.r.t. a selected reference configuration. The MCTDH[n] methods are able to treat large systems, but the single-reference Ansatz is not optimal in cases where one (or a few) degrees of freedom are special. Examples could be double-well systems, intramolecular vibrational-energy redistribution (IVR) calculations, or non-adiabatic dy-1 namics. In this work we introduce a multi-reference (MR) extension to the MCTDH [n] methods where selected higher-order excitations for the special degrees of freedom can be introduced in a simple but flexible way. The resulting MR-MCTDH[n] methods allow for e.g. treating non-adiabatic dynamics within the single-set formalism with the wave packets on each electronic surface described using the same level of approximation. Example calculations are performed on formyl fluoride (IVR), salicylaldimine (double well), and pyrazine (non-adiabatic dynamics). The results show that fast convergence is achieved by extending the configuration space in the special modes that govern the quantum dynamics.