The multiconfiguration time-dependent Hartree (MCTDH) method is formulated for density operators and applied to their numerical propagation. We introduce two types of MCTDH density operators which are expanded in different kinds of so-called single-particle density operators. The latter may either be hermitian, or else represent ket-bra products of so-called single-particle functions. For both types of MCTDH expansions of density operators we derive equations of motion employing the Dirac-Frenkel/MacLachlan variational principle. Further an alternative set of equations of motion for the second type of density operators is proposed, which is not based on a variational principle but derived by taking partial traces. We thus obtain three sensible approaches within the framework of the MCTDH method which differ in their performance and properties. We investigate these approaches and their properties analytically and numerically. Our numerical results refer to a model of vibronic-coupling dynamics in the pyrazine molecule representing coupled electronic states with four vibrational modes and two and three electronic states respectively. We analyze the closed-system dynamics for this model with temperature-dependent initial states. The influence of temperature on state populations, on correlation functions and on absorption spectra is discussed. We assess the numerical performance of two of the three approaches and find that both can be very efficiently applied to investigate the type of systems studied here.