The results obtained from the QM/MM simulations of the reaction mechanism of PET4 and PEF4 hydrolysis by PETase elucidate distinctive hallmarks that closely align with the experimental evidence, offering a finer mechanistic understanding with molecular detail. A first calculation of the FE profiles was performed at the DFTB3 level via WTMD simulations along the reaction coordinate obtained via HLDA. Our results showed that the acylation stage of PET4 proceeds through one single concerted step, crossing a FE barrier ≈ 21.0 kcal/mol high, and forming an acyl-enzyme intermediate with a reaction energy of ≈ 4.8 kcal/mol. The PET4 deacylation proceeds in a concerted manner as well, with activation and reaction energies of ≈ 12 kcal/mol and -1.7 kcal/mol, respectively, thus suggesting that the acylation is the rate-limiting step of the mechanism. These results agree both qualitatively and quantitatively with the umbrella sampling simulations performed at the PBE and DFTB3 level by Jerves et al. and García-Meseguer et al., respectively. However, in the work of Burgin et al., the FE barrier associated with the deacylation results ≈ 2 kcal/mol higher than that of the acylation, leading the authors to identify the acyl-enzyme breakdown as rate-limiting step of the process. A similar conclusion was reached by Dos Santos et al., who obtained a barrier of ≈ 20 kcal/mol for deacylation after a posteriori correction of semiempirical AM1 calculation at the M06-2X level. Globally, the magnitude of the FE barriers obtained in this work, as well as those shown by previous studies, are consistent with the kinetic measurements performed for BHET hydrolysis at 30 °C by PETase , and our data support the hypothesis that the acylation is the slowest stage of the reaction mechanism. In light of the agreement between the results obtained for PET4 with those previously reported with similar levels of theory, our approach was extended to the investigation of the reaction mechanism of PEF4 hydrolysis, for which detailed computational investigations are currently lacking in the literature. The FE curves obtained along the HLDA coordinates for the PEF4 hydrolysis are characterized by barriers of ≈ 10 kcal/mol and ≈ 5 kcal/mol for the acylation and deacylation steps, respectively. Therefore, akin to PET4 hydrolysis, the acylation remains the rate-limiting step of the process. On the other hand, at variance of the PET4 acylation, the formation of the PEF4 acyl-enzyme intermediate occurs through two steps: beyond the first FE barrier, a metastable tetrahedral intermediate is formed, which evolves towards the products after a second energy wall ≈ 2 kcal/mol high is overcome. This shape was confirmed by the PBE-level FE profiles calculated via MetaFEP and represents a crucial difference from PET4, indicating not only that the nucleophilic attack of S160 to the PEF4 ester group occurs faster than PET4, but also that a metastable tetrahedral intermediate is formed. Moreover, by applying MetaFEP, we pushed further the accuracy of the Hamiltonian employed to calculate the energetics accompanying the reaction using a wavefunction-based post-HF method (RI-MP2) and the popular DFT hybrid functional used for the description of organic molecules (B3LYP). Interestingly, these calculations produced almost identical FE curves, characterized by barriers of height compatible with those observed at the DFTB3 and PBE levels, although more negative energy differences were found for the products of both acylation and deacylation steps. Most importantly, even these high-accuracy calculations corroborated the mechanistic peculiarities of the PEF4 depolymerization obtained at the DFTB3 level, showing that the acylation proceeds with the formation of a metastable intermediate and this remains the rate-limiting step of the process. In addition, both the FE barriers and the reaction energy of the entire process obtained for PEF4 are lower than those obtained for PET4, indicating that the PEF depolymerization is both kinetically and thermodynamically favored. A possible explanation to this latter evidence is provided by the principles of conceptual DFT to compute key descriptors of the global and local reactivities of the two substrates. The global indices describing the electronegativity (𝜒), the chemical hardness (𝜂) and the electrophilicity (𝜔), suggested that PEF4 is more likely to accept additional electronic charge from the enzymatic environment. A deeper insight into in the local reactivities of the two substrates was given by the local electrophilicity (𝜔 ) and the local nucleophilicity (𝑁 6 ) defined from the condensed Fukui functions. These descriptors are sensitive to the method used to compute the atomic charges because of the inherent arbitrariness in density partitioning when defining an atom in a molecule. Moreover, a further factor that affects the estimation of the Fukui functions is the extension of the basis set. Indeed, the assessment of the local reactivities of the carbonyl group in different molecules conducted in the framework of the B3LYP theory showed some discrepancies among various atomic charge schemes dependent on the chosen basis set, with ChelpG being the only, notable exception. For these reasons, to mitigate the impact of any single method's limitations and provide a more robust understanding of reactivity patterns, it is advisable to select different partitioning schemes to conduct this analysis. From this standpoint, we compared the local reactivity indexes condensed over the atoms by means of the Fukui functions obtained via the Mulliken, ChelpG, Hirshfeld and MBIS methods, obtaining highly correlated patterns, descriptive of the different reactivity of the ester moieties in the two plastic chains. The increase of the 𝜔 6 function observed in PEF4 compared to PET4 indicate that ester group is more susceptible to the nucleophilic attack when located next to the furan ring rather than benzene. Conversely, the decrease of the 𝑁 6 9 function found for PEF4 suggests that this system is less reactive towards electrophilic reagents. This evidence explains the enhanced stability of the tetrahedral adduct formed upon S160 binding and its existence as a metastable intermediate. Indeed, we hypothesize that, in this state, the negative charge is more effectively delocalized in PEF4 compared to PET4. Furthermore, the ether O-atom (O5') in PEF4 is less prone to abstract the proton from the H237 histidine ring, which is the step that triggers the expulsion of the leaving group of the acylation and the formation of the acyl-enzyme intermediate. This characteristic is likely due to the electron-withdrawing effect of the O-atom of furan, that forms a dipole across the fivemembered ring and depletes the ester moieties of part of the electron density. In addition to the efficient catalysis, in our previous work we showed that PEF4 tends to dissociate quickly from the PETase binding pocket compared to its counterpart. The combination of the fast hydrolysis of PEF chains with the low residence time of the products of depolymerization points to a high turnover rate for the catalytic cycle, offering a valuable mechanistic description of the enhanced performance observed for PETase against PEF compared to its native substrate. Indeed, experimental findings obtained for macroscopic PEF and PET coupons depolymerization by PETase showed that, after 96 hours of incubation of the two polymers with the enzyme, the surface morphology of PEF resulted substantially more modified than that of PET. This observation was corroborated by differential scanning calorimetry (DSC) analysis, that showed a reduction in terms of relative crystallinity of 15.7 % for PEF compared with the relative reduction of 10.1 % for PET. Furthermore, recent analysis of Nano-PEF degradation by PETase showed that the WT enzyme is even more efficient than some popular variants engineered so far, which, in contrast, demonstrated highest catalytic activity towards PET. This finding can be ascribed to the stiffer active site that characterizes these mutants. While a rigid binding pocket enhances the thermostability of the enzyme, it also compromises the accommodation of the PEF chains, which are less flexible than those of PET because of the polarity of the furan ring and the nonlinear axis along which the aromatic ring is connected to the EG. On the other hand, other variants derived from PETase (FastPETase, TurboPETase 26 ) and the homologous LCC (LCC ICCG ) outperformed the WT enzyme in PET degradation and also demonstrated remarkable activity towards PEF. In line with the previous hypothesis, it is supposed that these enzymes retain a certain degree of the flexibility of the binding grooves, thus enabling remodeling of the catalytic domain to favor the accommodation of a broader range of substrates. Therefore, future engineering of these biocatalysts may be directed not only toward enhancing their performance in PET depolymerization, but also expanding their versatility towards other plastic materials.