  We present derivation and implementation of the Multi-Configurational
Strong-Field Approximation with Gaussian nuclear Wave Packets (MC-SFA-GWP) -- a
version of the molecular strong-field approximation which treats all electronic
and nuclear degrees of freedom, including their correlations,
quantum-mechanically. The technique allows, for the first time, realistic
simulation of high-harmonic emission in polyatomic molecules without invoking
reduced-dimensionality models for the nuclear motion or the electronic
structure. We use MC-SFA-GWP to model isotope effects in high-harmonics
generation (HHG) spectroscopy of methane. The HHG emission in this molecule
transiently involves strongly vibronically-coupled $^2F_2$ electronic state of
the $\rm CH_4^+$ cation. We show that the isotopic HHG ratio in methane
contains signatures of: a) field-free vibronic dynamics at the conical
intersection (CI); b) resonant features in the recombination cross-sections; c)
laser-driven bound-state dynamics; as well as d) the well-known short-time
Gaussian decay of the emission. We assign the intrinsic vibronic feature (a) to
a relatively long-lived ($\ge4$ fs) vibronic wave packet of the singly-excited
$\nu_4$ ($t_2$) and $\nu_2$ ($e$) vibrational modes, strongly coupled to the
components of the $^2F_2$ electronic state. We demonstrate that these physical
effects differ in their dependence on the wavelength, intensity, and duration
of the driving pulse, allowing them to be disentangled. We thus show that HHG
spectroscopy provides a versatile tool for exploring both conical intersections
and resonant features in photorecombination matrix elements in the regime not
easily accessible with other techniques.
