We present a theoretical framework for the calculation of rate constants of enzyme-catalyzed reactions that combines variational optimization of the dynamical bottleneck for overbarrier reactive events and multidimensional quantum mechanical tunneling dynamics for through-barrier reactive events, both in the presence of the protein environment. The theory features a two-zone, three-stage procedure called ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) with the transmission coefficient based on the equilibrium secondary-zone (ESZ) approximation for including the effects of the protein on a catalytic reaction center, called the primary zone. The dynamics is calculated by canonical variational theory with optimized multidimensional tunneling contributions, and the formalism allows for Boltzmann averaging over an ensemble of reactant and transition state conformations. In the first stage of the calculations, we assume that the generalized transition states can be well described by a single progress coordinate expressed in primary-zone internal coordinates; in subsequent steps, the transmission coefficient is averaged over a set of primary-zone reaction paths that depend on the protein configuration, and each reaction path has its own reaction coordinate and optimized tunneling path. We also present a simpler approximation to the transmission coefficient that is called the static secondary-zone (SSZ) approximation. We illustrate both versions of this method by carrying out calculations of the reaction rate constants and kinetic isotope effects for oxidation of benzyl alcoholate to benzaldehyde by horse liver alcohol dehydrogenase. The potential energy surface is modeled by a combined generalized hybrid orbital/quantum mechanical/ molecular mechanical/semiempirical valence bond (GHO-QM/MM/SEVB) method. The multidimensional tunneling calculations are inicrocanonically optimized by employing both the small-curvature tunneling approximation and version 4 of the large-curvature tunneling approximation. We find that the variation of the protein mean force as a function of reaction coordinate is quantitatively significant, but it does not change the qualitative conclusions for the present reaction. We obtain good agreement with experiment for both kinetic isotope effects and Swain-Schaad exponents.