QM regions included truncated heme-porphyrin, truncated substrates, and R242/245, H85/Q85 and two water molecules, the coordinates for the QM region can be found in SI. All protein residues and water molecules within 8Å of the Heme were included in the 'active region' of the QM/MM calculations. The atoms in the 'active region' interact with the QM atoms through electrostatic and van der Waals interactions and the corresponding polarization effects were considered in the subsequent QM/MM calculations. All QM/MM calculations were performed with ChemShell, by combining Turbomole for the QM part, and DL_POLY for the MM part. The MM region was described using the Amber ff14SB force field, and the electronic embedding scheme was used to account for the polarizing effect of the enzyme environment on the QM region. The QM/MM boundary was treated using hydrogen link atoms with the chargeshift model. During QM/MM geometry optimizations, the QM region was treated using the hybrid UB3LYP functional with two basis sets, B1 and B2, where B1 stands def2-SVP* and B2 is for def2-TZVP. For geometry optimization and frequency calculations, we used B1 basis set. The energy was further corrected using B3LYP-D3 37 functional. All of the QM/MM transition states (TSs) were located by relaxed potential energy surface (PES) scans followed by full TS optimizations using the P-RFO optimizer implemented in the HDLC code. The energies were further corrected with the large all-electron def2-TZVP basis set. The zero-point energies (ZPEs)