A scheme to efficiently calculate ro-vibrational (J > 0) eigenstates within the framework of the multiconfigurational time-dependent Hartree (MCTDH) approach is introduced. It employs a basis of MCTDH wave packets which is generated in the calculation of vibrational (J = 0) eigenstates via existing MCTDH-based iterative diagonalization approaches. The subsequent ro-vibrational calculations for total angular momenta J > 0 use direct products of these wave packets and the Wigner rotation matrices. In this ro-vibrational basis, the Hamiltonian matrix can be computed and diagonalized with minor numerical effort for any value of J. Accurate ro-vibrational states are obtained if the number of iterations in the J = 0 calculations and the basis set sizes in the MCTDH wave function representation are converged. Test calculations studying CH2D show that ro-vibrational eigenstates for moderately large J can be converged within wavenumber accuracy with the same MCTDH basis sets and only slightly increased iteration counts compared to purely vibrational (J = 0) calculations. If large J's are considered or very high accuracies are required, the number of iterations required to obtain convergence increases significantly. Comparing the theoretical results with experimental data for the out-of-plane bend, symmetric stretch, and antisymmetric stretch fundamentals, the accuracy of the ab initio potential energy surface employed is investigated.