Despite the remarkable achievements of QM/MM in practical applications, the majority of current approaches still treat nuclei in the key QM region classically, resulting in the neglect of nuclear quantum effects. This neglect is particularly problematic in systems where hydrogen, the lightest element, plays a significant role, as seen in many enzymatic reactions involving proton transfer, hydrogen atom transfer, and/or hydride transfer processes. To address this challenge, several theories have been developed to conduct QM/MM calculations with nuclear quantum effects included. One such successful approach is path-integral-based methods, which represent quantum systems using ensembles of replicas connected by harmonic springs. Pathintegral-based QM/MM methods have been successfully applied to study proton transfer, hydride transfer, and RNA cleavage reactions. Although a major limitation of path-integral based methods is the high computational cost, especially if ab initio potential energy surfaces for the QM part are to be used, recent developments have been able to accelerate the calculations and reduce their computational costs. Another promising approach is the nuclear-electronic orbital (NEO) method, which employs multicomponent wave functions to simultaneously describe the quantum behavior of both nuclei and electrons. When integrated with QM/MM, NEO has provided insights into the impact of nuclear quantum effects on molecular geometries in condensed phases. Although static NEO calculations are limited by the assumption of instantaneous quantum nuclear response to the motion of classical nuclei, the recently-developed real-time NEO QM/MM can address this limitation and incorporate nonadiabaticity between nuclei and electrons, thus offering insights into short-time vibronic dynamics. Additionally, semi-classical trajectory methods have been integrated with QM/MM calculations with nuclear quantum effects incorporated. They have been utilized to simulate the vibrational spectra of small biological molecules in condense phases. Recently, our group developed the constrained nuclear-electronic orbital (CNEO) theory to incorporate nuclear quantum effects, particularly quantum nuclear delocalization effects, into classical molecular simulations through a quantum-corrected effective potential energy surface. CNEO shows great potential to be a widely-used method for its simple physical picture, high computational efficiency, and accuracy for describing quantum nuclear delocalization effects. Due to its similarity to conventional electronic structure methods, with the addition of a more physically accurate quantum delocalized nuclear picture, CNEO is naturally capable of being integrated with the QM/MM framework.