Systemic engineering and global regulation enabling high-level bilirubin biosynthesis

Domain: physics
Published: 2025-12-23
Editorial outcome: major_revision

Abstract
Bilirubin biosynthesis has long been constrained by low yields and poorly understood bottlenecks. Here, we report a fully in vitro pathway that converts heme to bilirubin with the titer of 1.7 g/L and 95.8%. Systematically, enzyme screening and mechanistical analysis reveal the hidden challenge: Fe²⁺-induced oxidative degradation of intermediates. We show that Fe²⁺ coordinates with deprotonated intermediates to trigger oxidative ring-opening degradation via O₂-mediated radical mechanism, as supported by DFT calculations indicating reduced HOMO-LUMO gap in Fe²⁺-ligand complexes. The degradation is mitigated through competitive iron chelation and protonation state modulation, improving yield to 80.1%. Furthermore, we have resolved heme-CO complexes blocking O₂ activation at heme oxygenase by introducing a carbon monoxide dehydrogenase to remove CO and restore enzyme activity. Coupled with NADPH-recycling via formate dehydrogenase, these interventions enable efficient, scalable bilirubin synthesis with a 20-fold improvement. Our work shows controlling inhibitory byproducts is critical for stabilizing heme-related pathways and as a generalizable framework for synthetic biology.

Introduction
Bilirubin (BR), the principal bile pigment, was historically considered a metabolic waste product due to its pathological accumulation in conditions such as kernicterus and neonatal jaundice1–3. However, recent studies have revealed its beneficial physiological roles. BR possesses potent antioxidant activity, protecting lipids, proteins, and other biomolecules from oxidative stress4. Moderate serum BR levels are associated with reduced risks of cardiovascular and metabolic diseases3, and BR exerts neuroprotective effects by mitigating neuronal damage in the central nervous system5. Its anti-inflammatory, antioxidant, and antiviral properties have also been investigated in the context of COVID-19 treatment6. Furthermore, BR is the primary active compound in Calculus bovis, a traditional Chinese medicine with significant market value in China7.

Although BR is abundant in animal bile, conventional extraction from porcine or bovine sources is inefficient and cost-prohibitive, yielding only ~300 mg of BR per liter of bile8. Synthetic biology approaches have enabled enzymatic and de novo biosynthesis of tetrapyrrole compounds, including BR9,10. In nature, BR is derived from heme degradation: heme oxygenase (HO, EC 1.14.14.18) oxidizes heme to biliverdin (BV), which is subsequently reduced to BR by NADPH-dependent biliverdin reductase (BVR, EC 1.3.1.24)11–15 (Fig. 1a). HO isoforms exhibit regioselectivity by cleaving specific methylene bridges on the porphyrin ring to produce BV isomers (IXα–γ), but most BVRs selectively convert BV IXα into BR IXα—the biologically active form referred to hereafter as BR (Supplementary Fig. 1).

HO-mediated heme degradation requires three equivalents of O₂ and seven equivalents of NADPH, but HO lacks intrinsic cofactor-binding domains. Instead, electron transfer is facilitated by redox partner proteins (RPs), which typically consist of a large NADPH/FAD-binding domain and a smaller FMN-binding domain16–18. In some organisms, these domains are distributed across separate proteins and use Fe–S clusters to shuttle electrons19–21 (Fig. 1d). RP conformational changes regulate their activity: the closed conformation transfers electrons from NADPH to FMN, while the open conformation enables FMN to reduce the heme-bound HO22 (Fig. 1b). Subsequent O₂ binding and activation trigger porphyrin ring cleavage, generating BV, Fe²⁺, and carbon monoxide (CO) as byproducts (Fig. 1c, Supplementary Fig. 2)16.

BVR is a nicotinamide-dependent oxidoreductase with a Rossmann-fold N-terminal domain23,24. Structural studies have revealed that two BV molecules bind in the BVR active site; the distal BV may serve as an electron relay during the reduction of the proximal BV to BR24 (Fig. 1e and 1f; Supplementary Fig. 3). De novo heme biosynthesis via metabolic engineering has been widely explored25–28, and redirection of metabolic flux toward BR using HO and BVR has been attempted29. However, the complexity of the biosynthetic pathway has limited yields, with the highest reported BR titer at just 75.5 mg/L30. Alternatively, a two-step enzymatic conversion of heme to BR holds greater promise. Although diverse HO sources have been explored—including plant31, bacterial32, and cyanobacterial33 enzymes—BV titers remain modest (up to 132 mg/L)31, and BVR-mediated BR production has reached only 325 mg/L34. Notably, one-pot enzymatic synthesis of BR from heme has not yet been realized.

Importantly, the tetrapyrrolic structures of BV and BR retain metal-chelating properties similar to heme35. Transition metal complexes such as BV–Cu²⁺ are prone to oxidative degradation36,37, and BR can form Ca²⁺ and Cu²⁺ complexes that generate reactive oxygen species (ROS), leading to its breakdown38,39. In particular, Fe²⁺, a byproduct of HO activity, may compromise the structural integrity of both BV and BR. Additionally, CO—another byproduct—has strong ligand-binding properties. Its polarized bond and high electron density at the carbon center allow it to outcompete O₂ for coordination with heme iron40. CO binding stabilizes heme–CO complexes, which inhibit the formation of O₂-activated catalytic intermediates, thus preventing substrate access and further enzymatic turnover. These mechanistic insights suggest that uncontrolled byproduct interactions may be a critical factor limiting BR biosynthesis yields, which currently remain in the milligram per liter range.

In this work, we have constructed a compact enzyme library through phylogenetic analysis to identify high-activity variants of HO, RP, and BVR, which have been heterologously expressed in E. coli BL21 (DE3). Truncation of hydrophobic regions improves the soluble expression of HO and RP. To overcome the biosynthetic bottlenecks, we employ computational modeling to elucidate key mechanisms of byproduct interference. Density functional theory (DFT) simulations reveal Fe²⁺-induced O₂ activation and radical-driven oxidative degradation of BV/BR. On the other hand, we provide additional theoretical insights into the CO inhibition phenomenon, which is known in the community. CO-heme binding substantially increases the free energy of heme–CO complexes, inhibiting formation of O₂-coordinated intermediates and blocking substrate access to the enzyme active site. Guided by these insights, we develop an optimized bilirubin biosynthesis system that integrates metal chelation, protonation state control, and enzymatic CO elimination. This strategy enables efficient, gram-scale BR production, offering a robust platform for industrial-scale synthesis.

Results
Key enzyme selection for constructing the BR biosynthetic cascade
To establish an efficient in vitro BR biosynthetic pathway, we systematically selected and evaluated core enzymes responsible for the two-step conversion of heme to BR. In mammals and plants, two heme oxygenase (HO) isoforms—inducible HO-1 and constitutive HO-2—share high sequence and structural similarity, and a conserved catalytic mechanism41. Microorganisms typically encode a single HO isoform that catalyzes BR IXα formation, although some bacteria harbor HugZ-family heme-degrading enzymes that yield alternative isomers, including BR IXβ and BR IXδ42. Redox partner proteins (RPs) supporting HO activity also vary across taxa: animals and fungi employ NADPH-cytochrome P450 oxidoreductases (CPRs; EC 1.6.2.4)43, whereas plants and bacteria use a two-component system comprising either a ferredoxin (FDX) or flavodoxin as the primary electron acceptor and a corresponding ferredoxin/flavodoxin-NADP⁺ reductase (FDR) for electron donation20,21.

Biliverdin reductase (BVR), which catalyzes the dependent reduction of biliverdin (BV) to BR, is broadly distributed among animals and bacteria. In animals, BVRs fall into two classes—BVR-A and BVR-B subfamilies—distinguished by their regioselectivity44,45. Most Bacterial BVRs predominantly generate BR IXα and occasionally utilize F₄₂₀H₂ as a cofactor in place of NADPH46.

To identify the most effective enzyme combinations for the heme-BV-BR cascade, we conducted phylogenetic analysis of HO (Fig. 2a), RP (Fig. 2b), and BVR (Fig. 2c), selecting a panel of evolutionarily diverse candidates with potential for high catalytic efficiency. These enzymes were heterologously expressed in E. coli BL21 (DE3) (Supplementary Fig. 4), purified, and evaluated in vitro.

For the heme-to-BV conversion step, we tested various HO-RP pairings. Among HO-RP combinations, Rattus norvegicus HO-1 (RnHO-1) paired with its native CPR (RnCPR) exhibited the highest catalytic efficiency in converting heme to BV, producing 53.2 mg/L BV within 10 minutes (Fig. 2d). For the subsequent BV-to-BR conversion, the R. norvegicus BVR-A isoform (RnBVR-A, hereafter referred to as RnBVR) achieved the best performance, generating 338 mg/L BR in 10 minutes (Fig. 2e). These results established RnHO-1, RnCPR, and RnBVR as a high-performing triad for reconstituting BR biosynthesis in vitro.

Solubility-boosting truncations of HO and CPR retain catalytic performance
The expression of full-length of RnHO-1 and its redox partner RnCPR in E. coli was initially poor, consistent with the predictions of membrane-associating regions within their sequences. Structural analysis revealed that the N-terminal residues 1–54 of RnCPR and the C-terminal residues 268–289 of RnHO-1 as highly hydrophobic domains likely responsible for membrane anchoring (Fig. 3a), in line with prior studies implicating these regions in endoplasmic reticulum localization47. To enhance soluble expression, we implemented a progressive truncation strategy to remove these membrane-associated regions. Soluble expression levels increased in proportion to the truncation extent, with the fully truncated variants—RnHO-1Δ22 and RnCPRΔ54—exhibiting the highest expression levels (Fig. 3b).

We next evaluated whether truncation affected enzymatic performance. After optimizing the RnCPR-to-RnHO-1 ratio for maximal catalytic activity (Supplementary Fig. 5), we assessed the kinetic properties of the truncated variants. During testing, decreased catalytic rates of HO and its variants were observed at higher heme concentration range. Consequently, the Michaelis-Menten equation incorporating a substrate inhibition term was employed for fitting, yielding kinetic parameters including Km. All combinations of truncated RnHO-1 and RnCPR maintained comparable heme affinity (Km) and catalytic efficiency (kcat/Km) relative to their full-length counterparts (Fig. 3c and Supplementary Fig. 6), indicating that removal of membrane-associated domains did not compromise catalytic function. In contrast, strong substrate inhibition in the downstream reaction catalyzed by RnBVR prevented accurate kinetic modeling (Supplementary Fig. 7), underscoring the importance of enzyme-specific optimization strategies in cascade system design.

Byproduct as a hidden bottleneck in heme-based enzymatic cascades
Despite the use of cell lysates with equal concentrations of catalysts, the combination of RnHO-1Δ22 and RnCPRΔ54 exhibited a significantly faster reaction rate than the wild-type enzyme pair. However, an unexpected deviation from stoichiometric predictions was observed. Specially, while heme consumption and BV accumulation were anticipated to be directly proportional, at 80% substrate conversion, the BV yield only reached approximately 50%, with BV levels declining sharply as substrate consumption slowed (Fig. 4a). Notably, this issue was absent in the BVR-catalyzed BV-to-BR conversion, where reaction kinetics aligned with the expected pattern (Fig. 4b).

The faster reaction rate observed with RnHO-1Δ22/RnCPRΔ54 suggested high catalytic activity in the heme-to-BV step. However, the discrepancy in BV yield indicated the involvement of an additional factor influencing the reaction. Given that RnBVR demonstrated significantly higher specific activity (Fig. 3c), which catalyzed the subsequent BV-to-BR conversion, we hypothesized that the slower first step of the cascade (heme-to-BV) might be the bottleneck. To optimize reaction conditions and better balance cascade kinetics, we refined the experimental setup (Supplementary Fig. 8).

After optimizing conditions, we performed one-pot BR synthesis from heme. The heme conversion rate increased to 83.9% after 150 minutes, but the BR yield plateaued at 48.1% (Fig. 4c). This low yield could not be attributed to external factors, such as UV irradiation48 or ROS49 (e.g., H₂O₂), which are known to cause degradation of BR into byproducts like propentdyopents, as our reactions were conducted in the dark and without exogenous ROS. In contrast, the RnBVR-catalyzed reaction performed as expected, suggesting that the issue resides in the earlier heme-to-BV conversion step. We hypothesize that byproducts generated during this conversion, such as Fe²⁺ or CO, may destabilize BV or BR, thereby limiting the final BR yield.

Iron-induced activation of reactive intermediates drives oxidative degradation
We first explored the role of Fe²⁺, a byproduct of HO-catalyzed heme degradation, in promoting the reactivity and oxidative degradation of BR and its intermediate BV. Incubating a standard solution of BV and BR with Fe²⁺ under enzymatic reaction conditions (35°C, pH 7.5) for 150 minutes led to the emergence of multiple unknown compounds, as detected by HPLC (Fig. 5a and 5b). This was accompanied by a significant alteration in the UV-Vis spectra (Supplementary Fig. 9), which correlated with a reduction in the peak areas of BV and BR. LC-MS analysis identified 13 potential degradation products (Supplementary Figs. 10–12).

Given the established coordination chemistry between porphyrins and metal ions, we hypothesized that Fe²⁺ initially coordinates with BV/BR, thereby triggering degradation. DFT-based electrostatic potential surface analysis suggested that deprotonation of the central nitrogen atoms in BV and BR reduces their local electrostatic potential, enhancing Fe²⁺ coordination (Supplementary Fig. 13a). Time-dependent DFT (TD-DFT) calculations further indicated that deprotonation of BV induces redshifts in the Soret (~380 nm) and Q bands (690 nm), with an amplified intensity in the Q-band (Supplementary Fig. 13b and 13c). These changes aligned with the experimental UV-Vis alterations observed (Fig. 6c), confirming that Fe²⁺ preferentially interacts with deprotonated BV/BR (denoted as deBV/deBR), aligning with prior studies on porphyrin-metal ion coordination37.

DFT calculations revealed favorable coordination energies for deBV-Fe²⁺ (average –132.7 kcal/mol, Fig. 5c) and deBR-Fe²⁺ (average –153.4 kcal/mol, Fig. 5d). Notably, upon Fe²⁺ binding, deBR transitions from a saddle-shaped to a distorted planar conformation. LC-MS confirmed Fe²⁺-mediated oxidative degradation of BV and BR without the need for exogenous reactive oxygen species (ROS). We proposed that Fe²⁺ serves as a coordination core, recruiting axial ligands such as O₂ or OH⁻ from the solvent, thereby enhancing the reactivity of the entire complex. To explore this hypothesis further, we investigated various axial coordination modes (Supplementary Figs. 14 and 15). Among these, deBV-Fe-OH⁻ exhibited the lowest coordination energy. However, OH⁻ ligands poorly attack the porphyrin ring. Thus, plausible coordination modes, including deBV-Fe-O₂-OH (–11.8 kcal/mol) and deBR-Fe-O₂ (–26.1 kcal/mol), were identified. To prevent potential biases arising from a single functional/basis set combination, energy validation was performed across all spin states using three distinct functional/basis set combinations (Fig. 5c and 5d). Benchmark studies have demonstrated that different DFT methods exhibit distinct and systematic biases in predicting the spin states of transition metal systems50. Accordingly, the three selected functional/basis set combinations yielded different ground-state spin multiplicities across the systems studied. At all levels of theory, the ground state of deBV–Fe was consistently found to be a quintet. For deBV–Fe–O₂–OH, the ground state was a quintet with both the BP86-D3 and PBE0-D3 functionals, whereas it was a singlet when using DSD-PBEP86-D3. In the case of deBR–Fe, the quintet state was the ground state with BP86-D3, while the ground state was a triplet with both PBE0-D3 and DSD-PBEP86-D3. As for deBR–Fe–O₂, the ground state was a triplet with PBE0-D3, but a singlet with both BP86-D3 and DSD-PBEP86-D3. While PBE0-D3 calculations tend to overestimate the coordination free energies of Fe²⁺, all three DFT methods—based on their respective ground states—consistently indicate that coordination by deprotonated BV and BR, as well as subsequent axial ligand binding, is energetically favorable and potentially spontaneous.

Frontier molecular orbital (FMO) analysis revealed that Fe²⁺ coordination minimally alters the HOMO/LUMO distributions on the porphyrin ring (Fig. 5e and 5f). However, axial ligand binding (such as O₂) reduces the LUMO electron density on the porphyrin ring, shifting LUMO localization toward the O₂ ligand. Notably, Fe²⁺ coordination had minimal impact on the overall HOMO-LUMO gap, but the introduction of axial ligands significantly narrowed the gap (from 2.07 eV to 0.56 eV for deBV-Fe-O₂-OH and from 2.83 eV to 0.72 eV for deBR-Fe-O₂). These findings suggest that complexes with axial ligands, such as deBV-Fe-O₂-OH and deBR-Fe-O₂, exhibit enhanced reactivity, rendering them energetically favorable for single-electron transfer.

Given that FMO analysis indicated O₂ may act as an electron acceptor, and considering previous studies on heme enzyme-mediated ring-opening reaction mechanisms51, we proposed that this coordination facilitates the formation of superoxide radicals (O₂ → O₂⁻•), which subsequently trigger a radical cascade leading to meso-carbon attack and oxidative degradation (Fig. 5g). The proposed mechanism includes the following steps: (1) direct radical addition of O₂⁻• to meso-carbon, generating a Fe-O-O-carbon intermediate with a β-carbon radical; (2) homolytic cleavage of the O-O bond, forming epoxide intermediates and causing Fe-O bond cleavage; (3) hydrogen transfer from meso-carbon to the epoxide ring, resulting in the formation of a meso-carbon radical; and (4) free radical addition of another superoxide radical to meso-carbon, inducing a ring-opening reaction and leading to oxidative degradation products. This mechanism positions Fe²⁺ not only as a passive byproduct of heme degradation, but also as an active driver of porphyrin instability, promoting both structural disruption and radical initiation.

Preventing Fe²⁺-induced degradation via chelation and protonation control
To mitigate Fe²⁺-induced degradation of BV and BR, we systematically evaluated a panel of metal chelators. Although all tested chelators preserved the characteristics absorption of BV and BR (Supplementary Fig. 16), tartrate and other chelators compromised enzyme activity in both HO- and BVR-mediated reactions (Fig. 6a and 6b). Etidronic acid (HEDP) emerged as the most effective chelator, significantly enhancing the yields of BV and BR by 149.8% and 145.6%, respectively (Fig. 6a and 6b). Acid-base titration and DFT calculations confirmed that HEDP exists as di- and tri-deprotonated forms at physiological pH, and coordination titration established a 2:1 stoichiometry between HEDP and Fe²⁺ (Supplementary Fig. 17). Despite the highly favorable binding free energy (ΔG = –106.9 kcal/mol, Supplementary Fig. 18), residual degradation products will still detectable (Fig. 6e and 6f), suggesting that the suppression of deBV/deBR-Fe²⁺ complex formation was incomplete.

Given that Fe²⁺ coordination is strongly pH-dependent, we next investigated whether protonation control could enhance suppression of degradation. While BR readily precipitated under acidic conditions, BV remained soluble and exhibited progressive redshifts in the Soret and Q bands with increasing pH (3.5‒9.5), consistent with stepwise deprotonation (Fig. 6c). When Fe²⁺ was added under alkaline conditions drastically, BV absorption dropped sharply, whereas its spectral profile was largely preserved under acidic conditions (Fig. 6d), suggesting that lower pH impedes Fe²⁺ coordination to deprotonated BV. Combining HEDP chelation with mild acidification (pH 6.5), resulted in a synergistic reduction of BV and BR degradation (Fig. 6e and 6f). This optimized condition, including 15 mM HEDP and pH 6.5, preserved enzyme activity while minimizing Fe²⁺-induced side reactions and maximizing product yields (Supplementary Fig. 19). These results underscored the importance of both chelation and pH tuning in maintaining the stability of heme-derived products under enzymatic conditions.

Effects of gaseous byproduct: verification of CO inhibition
Under optimized HEDP dosage and pH conditions to suppress Fe²⁺ interference, the BR yield increased significantly to 80.1% (1953 mg/L, 3 mM heme as substrate) with the RnHO-1Δ22, RnCPRΔ54, and RnBVR mixture. Despite this enhancement, ~20.5% of the heme remained unconverted during the reaction’s latest stages (Fig. 7a), suggesting the presence of an additional bottleneck. Given the known inhibitory properties of CO—a stoichiometric byproduct of heme oxygenase (HO) catalysis—assessment of CO inhibitory effects on the established heme oxygenase system is required52. CO is generated from the meso-carbon atom of heme and released via a conserved gas channel within the HO enzyme (Fig. 7c).

Employing three distinct levels of theory, including the double-hybrid functional DSD-PBEP86-D3, it was confirmed that CO forms a thermodynamically more stable complex with heme compared to O₂ (Fig. 7d). This supports CO’s role as a potent competitive agent capable of displacing O₂ from its coordination site. Although the spin of the ground state differed under various computational conditions, the average binding energy of CO to Heme (–23.2 kcal/mol) was significantly stronger than that of O₂ (–6.9 kcal/mol). Supporting this notion, CO-pretreated heme solutions showed markedly diminished HO activity (Fig. 7b), indicating effective inhibition of catalysis.

To further verify the molecular basis of CO inhibition on enzyme activity, we performed molecular dynamics simulations of HO complexes with heme-O₂ and heme-CO. The computed binding free energies revealed a stronger enzyme affinity of heme-O₂ (–2.50 kcal/mol) than for heme-CO (5.54 kcal/mol) (Fig. 7e). Per-residue decomposition analysis of the binding free energies identified key residues around the substrate-binding pocket that contribute to this difference (Fig. 7f and 7g). Notably, Ser142, located adjacent to the CO ligand-binding site, exhibited higher steric energy levels in the heme-CO complex compared to heme-O₂. This arises from the rigid, linear geometry of the sp-hybridized CO ligand, which created spatial clashes in contrast to the more flexible, bent sp²-hybridized O₂ ligand. This geometric disparity induces steric hindrance consistent with crystallographic observations53, ultimately resulting in significantly lower affinity of HO for heme-CO than for heme-O₂.

The energy contribution comparison between heme-CO and heme-O₂ complexes within the substrate pocket residues (Fig. 7h) highlighted widespread energetic penalties in the CO-bound state. Furthermore, energy landscape mapping based on RMSD and radius of gyration during stable simulation phases revealed a deeper, more stable conformational basin for the heme-O₂ complex (Fig. 7i, Supplementary Fig. 20), whereas the heme-CO complex was associated with a shallower and less favorable energy minimum (Fig. 7j).

Collectively, these results provide a supplementary verification of structural and thermodynamic rationale for CO-mediated inhibition: CO acts as a stronger ligand than O₂, forming a stable heme–CO complex with free heme in the solution. However, pre-formed complex exhibits steric incompatibility with the HO catalytic pocket. Even if accommodated, CO occupies the axial position of the iron center, where its strong coordination prevents O₂ displacement and activation, thereby halting the subsequent catalytic cycle.

Integrated enzymatic cascade enables gram-scale bilirubin production
To overcome CO-mediated inhibition of HO, we incorporated an engineered carbon monoxide dehydrogenase (ChCODH)54 to catalyze the oxidization of CO into non-toxic CO₂. Concurrently, an NADPH-dependent formate dehydrogenase (ApFDH)55 was introduced to enable in situ cofactor regeneration, the reducing equivalent essential for the HO and BVR reactions. Therefore, building upon the reaction progress depicted in Fig. 7a, we evaluated the BR yield by supplementing the system with cell lysate containing either ApFDH or ChCODH. The results demonstrated that supplementation with ApFDH achieved the same yield as stoichiometric NADPH addition, without the need for external NADPH. Moreover, when ChCODH was added, the effective removal of CO led to a further increase in the BR yield from 80.3% to 91.3% (Fig. 8a). For efficient multi-enzyme co-expression, we employed pCDFDuet-1 and pETDuet-1 plasmids, supplementing the latter with an additional ribosome binding site (RBS) to balance translation efficiency. This assembly yielded a five-enzyme co-expression system comprising RnHO-1Δ22, RnCPRΔ54, RnBVR, ApFDH, and ChCODH (Supplementary Fig. 21), forming a fully integrated multi-enzyme catalysis cascade (Fig. 8b).

Using this system, we achieved complete conversion of 3 mM (1953 mg/L) heme to 2.87 mM (1678.4 mg/L) BR, corresponding to a 95.8% yield in a 1 L in vitro reaction (Fig. 8d). The reaction was visually tracked by a color change from dark brown to reddish-brown (Fig. 8c). Following extraction and purification using rotary evaporation and lyophilization method, the powder product was confirmed as the BR-IXα isomer via 1H and 13C NMR spectroscopy (Supplementary Fig. 22 and 23).

This study represented a gram-scale biosynthesis of BR through a fully reconstituted, enzyme-based system (Fig. 8e). While existing research on bilirubin biosynthesis often utilizes readily available carbon sources such as glucose, which offers cost-effectiveness, our cascade system achieved a 20-fold higher yield compared to prior microbial or enzymatic systems30. Unlike systems relying on carbon feedstocks such as glucose, our optimized in vitro platform bypassed metabolic burden while achieving high efficiency. More importantly, our end-to-end cost analysis for the developed process indicates that BR can be produced at 0.22 g per dollar, beginning from cell cultivation. This cost efficiency is highly competitive with de novo synthesis approaches. Furthermore, the rational truncation of membrane-associated domains to enhance soluble protein expression, coupled with the systematic elimination of inhibitory byproducts, provided a modular framework for engineering redox-intensive biosynthetic pathways.

Methods
Materials
All wild-type gene sequences were codon-optimized, synthesized, and cloned into the BamHI/XhoI restriction sites of the pET28a(+) vector by Sangon Biotech (Shanghai, China). PrimeSTAR Max DNA Polymerase and the restriction enzyme DpnI were obtained from Takara Bio Inc. (Japan). The homologous recombination and DNA purification kits were purchased from Vazyme Biotech Co., Ltd. (Shanghai, China). The SDS-PAGE gel preparation kit was sourced from Sangon Biotech (Shanghai, China). NADPH was acquired from Roche Holding AG (Switzerland), and heme (hemin chloride, molecular weight = 651.94 Da) was purchased from Aladdin Industrial Corporation (Shanghai, China). Plasmids pETDuet-1, pCDFDuet-1, and E. coli BL21 (DE3) were preserved in our laboratory stock. All other reagents were sourced from Shanghai Macklin Biochemical Co., Ltd. (Shanghai, China), unless otherwise specified.

Phylogenetic analysis
Phylogenetic trees were constructed using the Neighbor–Joining method in Molecular Evolutionary Genetics Analysis version 11 (MEGA11)59. Protein sequences used for tree construction were retrieved from the NCBI database, with selected sequences (Supplementary Table 1) representing diverse biological origins and catalytic properties. The tree was visualized using the Chiplot online platform60.

Protein expression and purification
Recombinant plasmids were transformed into E. coli BL21 (DE3) strains for protein expression. E. coli cultures were grown in LB medium (10 g/L tryptone, 5 g/L yeast extract, 5 g/L NaCl) with kanamycin or ampicillin (50 μg/mL) at 37°C with shaking (200 rpm). Protein expression was induced by adding isopropyl-β-D-thiogalactopyranoside (IPTG) to a final concentration of 0.1 mM when the OD600 reached 0.4‒0.6. After 14 hours of induction at 20°C, cells were collected by centrifugation (8000 × g), and SDS-PAGE analysis was performed to verify protein expression levels. The enzyme was further purified using an ÄKTA Protein Purifier system, cells were resuspended in 100 mM Tris-HCl buffer (pH 7.5), and disrupted by ultrasonication on ice. The supernatant was obtained by centrifugation (12,000 × g for 10 min at 4°C) and loaded on a HisTrap HP chromatographic column (GE Healthcare, Chicago, USA) with buffer A (100 mM Tris-HCl, 150 mM NaCl, pH 7.5), and the enzyme was eluted with buffer B (100 mM Tris-HCl, 150 mM NaCl, 0.5 M imidazole, pH 7.5). The target protein solution was collected and concentrated with an ultrafiltration tube (Millipore Amicon Ultra, Billerica, USA), and then a desalting column (Merck KGaA, Darmstadt, Germany) was used to remove imidazole61.

Enzyme screening
Reactions for HO and its partner proteins, CPR or FDR+FDX combination, were conducted by adding 100 mg/L HO and either 100 mg/L CPR or 100 mg/L FDR + 100 mg/L FDX, with 1 mM heme and 10 mM NADPH, followed by incubation at 35°C with agitation (600 rpm) for 10 min. Reactions were quenched by a 10-fold dilution with DMSO, and product concentrations were analyzed by HPLC. For BVR screening, reactions containing 100 mg/L BVR, 1 mM BV, and 2 mM NADPH were conducted at 35°C with agitation (600 rpm) for 10 min, followed by quenching with DMSO dilution, and product levels were determined by HPLC.

Generation of truncated mutants
Truncated plasmids were generated by amplifying gene sequences with the primer pairs listed in Supplementary Table 2, using pET28a(+)-RnHO-1 or pET28a(+)-RnCPR as templates, with the following PCR protocol: the PCR reaction system was carried out in 50 μL of a solution consisting of 25 μL DNA PrimeSTAR Max DNA polymerase, 1 μL F- and R-primer, 1 μL DNA template, and 22 μL ddH2O. The PCR program was 98°C for 30 s, 55°C for 1.5 min, 72°C for 10 min, a total of 30 cycles62. The PCR products were treated with DpnI to remove template plasmids, purified, and subjected to homologous recombination using a Vazyme Biotech kit. The resulting plasmids were transformed into E. coli BL21 (DE3), and colonies were screened for the correct mutants, which were validated by DNA sequencing (Sangon Biotech, Shanghai, China).

LC-MS analysis
LC-MS analysis of BV and BR degradation products was performed using a Dionex UltiMate 3000 UHPLC system (Thermo Fisher Scientific, Waltham, MA, USA) coupled with a Thermo TSQ Quantum Ultra triple quadrupole mass spectrometer. Chromatographic separation was achieved using a Waters ACQUITY UPLC BEH C18 column (100 mm × 2.1 mm, 1.7 μm particle size; Waters Corporation, Milford, MA, USA) with electrospray ionization (ESI) in positive ion mode. Mobile phase A, water (0.1% formic acid); mobile phase B, acetonitrile (0.1% formic acid). The gradient elution program was as follows: 0‒5 min: 20% B; 5‒10 min: 20‒40% B; 10‒15 min: 40% B; 15‒20 min: 40‒100% B; 20‒30 min: 100% B; 30‒32.5 min: 100‒20% B; 32.5‒35 min: 20% B. The column temperature was maintained at 30°C, and the flow rate was 0.2 mL/min. Data were processed using Xcalibur software (Thermo Fisher Scientific).

Screening of chelators
For initial chelator screening, 1 mM BV or BR standard solutions were mixed with 5 mM chelator candidates (tartrate, citrate, EDTA, DEG, HEDP, NTA, HEDTA) and 1 mM FeCl₂. After incubation at 35°C for 150 minutes, samples were analyzed by UV-Vis spectroscopy (Supplementary Methods). To evaluate chelator effects on enzymatic activity, crude enzyme reactions were conducted with RnHO-1Δ22, RnCPRΔ54, or RnHO-1Δ22, RnCPRΔ54, and RnBVR in the presence of 5 mM of each chelator. Reactions were performed under conditions (1 mM heme, 10 mM NADPH, 100 mM Tris-HCl, pH 7.5, 35°C, 600 rpm), and products were quantified by HPLC.

Theoretical calculations
Density functional theory (DFT) calculations were performed using Gaussian 09, Revision D.0163. Initial molecular structures were modeled in GaussView 5.0 and subjected to conformational searches with the Molclus program64. Geometry optimizations and vibrational frequency analyses were performed at the B3LYP/6-311G(d,p) level of theory, and single-point energy calculations were carried out at the BP86/def2-TZVP level. For key species, calculations were additionally performed using two DFT methods: PBE0/def2-TZVP and DSD-PBEP86/def2-TZVPP. A polarizable continuum model (PCM) was used to simulate aqueous solvation effects, and Grimme’s GD3BJ dispersion correction was applied. For Fe²⁺-containing systems, all possible spin multiplicities were evaluated to determine the ground-state configuration. UV-Vis absorption spectra for BV were calculated using time-dependent DFT (TD-DFT) at the BP86/def2-TZVP level, with 50 excited states computed. Electrostatic potential (ESP) charges and frontier molecular orbitals (FMOs) were analyzed using the Multiwfn program (Version 3.8)65, and results were visualized using VMD software (Version 1.9.3)66.

Scaled-up BR synthesis
For scaled-up BR synthesis, the reaction was conducted in a 5 L bioreactor (T&J Bioengineering, Shanghai, China). The system contained 20 g dry cell weight (DCW)/L of lysate from the five-enzyme co-expression strain E. coli/pCDFDuet-1-RnCPRΔ54/RnHO-1Δ22/pETDuet-1-ChCODH/RnBVR/ApFDH, 1953 mg/L (3 mM) heme, 50 mM ammonium formate, 0.1 mM FMN, and 15 mM HEDP. The reaction was conducted at 35°C with agitation (400 rpm), with pH automatically maintained at 6.5 ± 0.1 by controlled addition of HCl, and air was continuously supplied at 2 vvm (volume per volume per minute).

Cell lysate for the 1 L scale was prepared using a high-pressure homogenizer (ATSAH-1500) operated at 850 bar with a flow rate of 5 L/h. The cell suspension was passed through the homogenizer three times to obtain the lysate, with chilled circulating water maintained at 4°C.

Statistical analysis
Each experiment was performed at least in triplicate. Values were expressed as mean ± standard deviation (SD). Statistical analysis was performed by using Origin 2024 with a two-tailed unpaired Student’s t-test. And a value of p < 0.05 was considered statistically significant.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present a one pot in vitro enzyme cascade converting heme to bilirubin (BR) with g/L titers and ~95% yield, surpassing the previous best published pathways by an order of magnitude (based on volume of cell growth versus volume of cell-free reaction). This work systematically screens across enzymes and utilizes computation guided mechanistic analysis to understand reaction dynamics, predict and address known and unknown bottlenecks, and propose enzyme mechanisms that are supported by collected data. There are 2 major bottlenecks addressed and overcome in this paper which include iron induced degradation of reaction intermediates, and inhibited heme oxygenase activity by carbon monoxide. The authors do a great job of addressing each step in the pathway and associated bottlenecks with computationally supported hypothesis and clear mechanistic insights. Experimental results confirm the production of BR at the gram scale in a 1L bioreactor. There are a few remaining questions about some of the work but overall, the authors present a compelling story displaying how fundamental chemical insights can guide pathway engineering for production of an important chemical product at the g/L scale. While a good story, I believe the authors need to tone down some of the verbiage associated to \"record\", \"transcends\" , etc. Let the data speak for themselves. In regards to the data, better comparisons to state of the art need to be carried out as I describe in my second major point below.\nMajor comments:\n1. The authors provide significant experimental and computational data to suggest the CO mediated inhibition of heme oxidase and use MD simulations to identify key residues in the substrate binding pocket that sterically clash with the heme-CO complex. To address this claim, the authors add a CO dehydrogenase enzyme and concurrently add an NADPH dependent formate dehydrogenase for cofactor regeneration. These additions improve yield and titer; however, it is not clear how these 2 variables (CO dehydrogenase and formate dehydrogenase) contribute individually to the observed yield and titer improvements. To confirm the author’s claim that CO is truly inhibitory, the authors need to test these two additions separately as controls at a small scale before implementing in 1L bioreactors.\n2. While the authors compare volumetric yields, additional comparisons are needed. To more convincingly show their approach is better than the state of the art, the authors need to do two things. (i) Please compare product yield relative to g cell pellet. In other words, how much cell pellet was needed to get the enzymes used in the in vitro cascade. This would be a better comparison to the in vivo results. The other way to think of it is to put everything on a volume basis of cells grown (i.e., cells grown to get to the enzymes needed for the in vitro cascade). (ii) The authors should present a simple cost model of end to end costs starting with cell growth for the yields achieved (e.g., g product per $). It may be that the in vitro cascade has a higher yield in the in vitro reaction, but is overall more costly, or is more expensive in terms of how many cells are needed to get to the final result. This needs to be calculated and discussed.\nMinor comments:\n1. Figure 1A: It is not clear why other BV intermediates are shown here, as it distracts from the main text. We suggest that these intermediates go into supplemental and the main 95% abundant intermediate is shown.\n2. SI fig 3: It is not clear why expression profiles are shown here as opposed to purified proteins, as all of these homologs were purified and tested in vitro.\n3. SI fig 5 and 6: Is it fair to disregard the last 3 points in each graph? Can you either: show what forecasted kinetics would look like without substrate inhibition (curve fitted) beyond 0.4 mM substrate or use Michaelis-Menton kinetic equations with substate inhibition terms to see how these compare\n4. Fig2A-C: these figures are not directly cited in the manuscript.\n5. Line 322: please delete the duplicated word (including)\n6. In the methods, could you provide how the lysate was made for the bioreactor?", "reviewer_2": "The study sits at the intersection of synthetic biology, biocatalysis and enzyme engineering, where the authors reconstitute or redesign multi-enzyme pathways to manufacture high-value natural products (the bile pigment bilirubin) from inexpensive precursors. By combining five enzymes in a single, cell-free condition, the authors increase bilirubin (a bile-acid precursor) titer from the long-standing ceiling of tens of milligrams per litre to 1.68 g L⁻¹ in <3 h. This is a >20-fold jump that turns a niche diagnostic dye into a scalable valuable product. This breakthrough demonstrates how pathway conditioning can unlock high-value tetrapyrrole manufacturing in synthetic biology. The study couples physics-based modeling (DFT + MD) with experimentally-driven cascade design. Quantum chemistry calculations detail how two coproducts of heme oxygenase —Fe²⁺ and CO—trigger biliverdin/bilirubin degradation or enzyme stalling; parallel bench work implements iron chelation, CO scavenging, and NADPH recycling to avert those bottlenecks, delivering gram-scale output. In general, I consider the study significant and innovative. Here are some comments:\n1. Relevant to the current study, a recent perspective focuses on the impact of physics-based modeling (QM and MD) on enzyme engineering. The present work illustrates a complementary aspect—intermediate-guided engineering—showing that managing unstable or inhibitory pathway intermediates (Fe²⁺, CO) can be more decisive for yield than enzyme optimization. I would suggest the author to conduct a super brief comparison to the perspective by Jurich et al. (“Physics-based modeling in the new era of enzyme engineering,” Nat. Comput. Sci. 5, 279–291, 2025), highlighting this contribution as a complementary design tier for synthetic biology.\n2. Major conceptual critique: computation is confirmatory, not hypothesis-generating. Before this study the community already knew that 1) free Fe²⁺ catalyses oxidative cleavage of biliverdin/bilirubin, and chelators or acidification can suppress the damage; 2) CO binds heme far more tightly than O₂, forming a photolabile heme-CO complex that inhibits heme-oxygenase turnover; employing CO-scavenging enzymes is a recognised workaround. It is true that the authors’ DFT and MD simulations further quantify these effects, but I would argue that these computations are primarily confirmatory; the chelator (HEDP) and CODH choices appeared to be guided by prior empirical knowledge rather than by de-novo quantum chemistry insight. I would suggest the authors to: 1) emphasize the design principle that the DFT results （homo-lumo gap and binding free energy） are uniquely poised to inform (some guidance by DFT that can not be obtained from simple empirical understanding); and/or 2) consider framing DFT more as a supportive approach to achieve the goal, rather than guiding the development of the approach.\n3. Other technical details. Spin-state coverage: only the quintet state is reported. It's important to consider the triplet state as it may lie near-degenerate for Fe²⁺ porphyrins.\n4. Functional/basis set validation. Functional/basis set – B3LYP-D3(BJ)/6-31G(d,p) is used here. Suggest to consider other range-separated functional and also DLPNO-CCSD(T) single-points on key species. This would enhance reliability.\n5. Limited sampling from MD. A single 200 ns replica of HO-1 is short for a membrane-associated protein; multiple replicas (3 x 200 ns) or enhanced sampling would firm up the Ser142 steric-block hypothesis.\n6. Data availability. In addition to the Cartesian coordinates, sharing input decks, force-field parameters, and MD trajectories will enhance the reproducibility of the current work."}, "author_response": "Response to reviewers’ comments\nReviewer #1:\nThe authors present a one pot in vitro enzyme cascade converting heme to bilirubin\n(BR) with g/L titers and ~95% yield, surpassing the previous best published pathways\nby an order of magnitude (based on volume of cell growth versus volume\nof cell-free reaction). This work systematically screens across enzymes and utilizes computation\nguided mechanistic analysis to understand reaction dynamics, predict and address\nknown and unknown bottlenecks, and propose enzyme mechanisms that are supported\nby collected data. There are 2 major bottlenecks addressed and overcome in this paper\nwhich include iron induced degradation of reaction intermediates, and inhibited heme\noxygenase activity by carbon monoxide. The authors do a great job of addressing each\nstep in the pathway and associated bottlenecks with computationally supported\nhypothesis and clear mechanistic insights. Experimental results confirm the production\nof BR at the gram scale in a 1L bioreactor. There are a few remaining questions about\nsome of the work but overall, the authors present a compelling story displaying how\nfundamental chemical insights can guide pathway engineering for production of an\nimportant chemical product at the g/L scale. While a good story, I believe the authors\nneed to tone down some of the verbiage associated to \"record\", \"transcends\", etc.\nLet the data speak for themselves. In regards to the data, better comparisons to state of\nthe art need to be carried out as I describe in my second major point below.\nResponse: We sincerely thank the Reviewer for this comment. We have moderated\nthe language throughout the manuscript, particularly terms such as \"record\",\n\"transcends\", and similar phrases, to ensure the clarity and an appropriate tone.\nThe specific revision are as follows:\nLine 99: “This strategy enabled gram-scale BR production with the highest yield\nreported to date” → “This strategy enabled efficient gram-scale BR production”;\nLine 440: Deleted “the highest reported titer to date”;\nLine 448: Deleted “remarkable”;\nLine 471: Deleted “record-breaking”;\nLine 520: “transcends incremental advances” → “provides a significant advance”\nwas deleted;\nThe entire manuscript has been rechecked to unify the tense as well as revise\nseveral formatting details, all of which have been marked in red.\nMajor comments:\n1. The authors provide significant experimental and computational data to suggest the\nCO mediated inhibition of heme oxidase and use MD simulations to identify key\nresidues in the substrate binding pocket that sterically clash with the heme-CO complex.\nTo address this claim, the authors add a CO dehydrogenase enzyme and concurrently\nadd an NADPH dependent formate dehydrogenase for cofactor regeneration. These\nadditions improve yield and titer; however, it is not clear how these 2 variables (CO\ndehydrogenase and formate dehydrogenase) contribute individually to the observed\nyield and titer improvements. To confirm the author’s claim that CO is truly inhibitory,\nthe authors need to test these two additions separately as controls at a small scale before\nimplementing in 1L bioreactors.\nResponse: Thank Review very much! In response to your suggestions, we have\ncomplemented the experiments at the 10 mL scale to validate the effects of the ApFDH\nand ChCODH enzymes. Specifically, the addition of ApFDH eliminated the need for\nexogenous NADPH supplementation, significantly reducing production costs. The\ninclusion of ChCODH increased the final bilirubin yield by approximately 11%,\nreaching 91.3%. These results further clarified the specific roles of these auxiliary\nenzymes within the main cascade pathway.\nWe have revised manuscript as following:\n(i) Fig. 8: Now includes the cascade reaction results with the addition of ApFDH\nand ChCODH. The figure layout has been adjusted, and the figure caption has been\nrevised according to reflect these changes. The revised figure was marked in yellow.\n2\n(ii) We have added the following description “Therefore, building upon the\nreaction progress depicted in Fig. 7a, we evaluated the BR yield by supplementing the\nsystem with cell lysate containing either ApFDH or ChCODH. Supplementation with\nApFDH achieved the same yield as stoichiometric NADPH addition, eliminating the\nneed for external NADPH. Moreover, the addition of ChCODH, which effectively\nremoved CO, increased the BR yield from 80.3% to 91.3% (Fig. 8a)”. (Line 426-432)\nWe have supplemented the following method description in Supporting\ninformation “Effects of ApFDH and ChCODH on the conversion of heme to BR.\nThe conversion of heme to BR was primarily catalyzed using cell lysates (each 20 g\ndry cell /L) of RnHO-1 , RnCPR and RnBVR. Reactions were performed in a 10\nΔ22 Δ54\nmL system containing 3 mM heme, 30 mM NADPH, 0.1 mM FMN, 15 mM HEDP,\nand 100 mM Tris-HCl (pH 6.5) at 35°C with 600 rpm agitation for 150 min.\nTo investigate specific effects on the reaction:\nNADPH recycling effect: 20 g/L cell lysate of ApFDH was supplemented while\nomitting the initial 30 mM NADPH addition.\nCO elimination effect: 20 g/L cell lysate of ChCODH was supplemented.\nSynergistic effect of ApFDH/ChCODH: Both ApFDH and ChCODH crude lysates\n(each 20 g/L) were co-supplemented, concurrently omitting the 30 mM NADPH\n3\naddition.” (Line 131-142).\nAdditionally, we believe that the data in Fig. 7b—specifically, the substantially\ndecrease in catalytic efficiency observed when using CO-saturated heme as the\nsubstrate—provide compelling evidence that CO significantly inhibits the heme\noxygenase (HO)-mediated oxidative ring-opening reaction of heme.\n2. While the authors compare volumetric yields, additional comparisons are needed. To\nmore convincingly show their approach is better than the state of the art, the authors\nneed to do two things. (i) Please compare product yield relative to g cell pellet. In other\nwords, how much cell pellet was needed to get the enzymes used in the in vitro cascade.\nThis would be a better comparison to the in vivo results. The other way to think of it is\nto put everything on a volume basis of cells grown (i.e., cells grown to get to the\nenzymes needed for the in vitro cascade). (ii) The authors should present a simple cost\nmodel of end to end costs starting with cell growth for the yields achieved (e.g., g\nproduct per $). It may be that the in vitro cascade has a higher yield in the in vitro\nreaction, but is overall more costly, or is more expensive in terms of how many cells\nare needed to get to the final result. This needs to be calculated and discussed.\nResponse: Thank Review very much for the thoughtful and constructive\nsuggestions! We have made point-to-point revision as follows:\n(i) Fig. 8e: The figure now includes the specific yield—defined as the mass of\nbilirubin obtained per unit cell dry weight (mg/g dry cells) —calculated for all studies\nreported in the literature. These results demonstrate that our method not only achieves\nthe highest reported yield to date, but also delivers the highest bilirubin titer under\nequivalent biomass conditions.\n(ii) End-to-end cost analysis: We calculated the total production costs from cell\ncultivation through the completion of the cascade reaction, expressed as grams of\nbilirubin per dollar ($.) This result (0.22 g bilirubin per dollar) has been added as:\n“More importantly, our end-to-end cost analysis for the developed process\nindicates that BR can be produced at 0.22 g per dollar, staring from cell cultivation.\n451-454)\nMinor comments:\n1. Figure 1A: It is not clear why other BV intermediates are shown here, as it distracts\nfrom the main text. We suggest that these intermediates go into supplemental and the\nmain 95% abundant intermediate is shown.\nResponse: We sincerely appreciate the reviewer's valuable suggestions! We have\nimplemented the recommended changes and updated Fig. 1a (Line 101) in the\nmanuscript. The original version of Fig. 1a has been moved to the Supplementary\ninformation, renamed as Supplementary Fig. 1, and cited in the text “but most BVRs\nselectively convert BV IXα into BR IXα—the biologically active form referred to\nhereafter as BR (Supplementary Fig. 1).”. (Line 51)\n2. SI fig 3: It is not clear why expression profiles are shown here as opposed to purified\nproteins, as all of these homologs were purified and tested in vitro.\nResponse: Thank Reviewer very much! According to Review’s good suggestions,\nwe have updated Supplementary Fig. 3 (now renumbered as Supplementary Fig. 4)\nto display electrophoresis bands of all homologs purified by Ni2+ affinity\nchromatography. The following description to the figure caption:\n“a, heme oxygenases; b, redox partners; c, biliverdin reductases. For samples\nwith poor solubility, Amicon® Ultra-15 centrifugal filter (10000 NMWL, Merck\nMillipore Ltd.) was used to concentrate them before conducting activity assays”.\n5\n3. SI fig 5 and 6: Is it fair to disregard the last 3 points in each graph? Can you either:\nshow what forecasted kinetics would look like without substrate inhibition (curve fitted)\nbeyond 0.4 mM substrate or use Michaelis-Menton kinetic equations with substate\ninhibition terms to see how these compare.\nResponse: We thank the Reviewer for the constructive suggestions! In the original\nSupplementary information (SI) Fig. 5, we incorrectly applied the conventional\nMichaelis-Menten equation without accounting for substrate inhibition at higher\nconcentrations. In the revised manuscript, we have:\n(1) Reanalyzed the data using the substrate inhibition-modified Michaelis-Menten\nequation:\n𝑣\n𝑚𝑎𝑥\n𝑣 =\n𝐾 [𝑆]\n1+ 𝑚 +\n[𝑆] 𝐾\n𝑖\n(2) Conducted additional independent replicates for datasets with large standard\ndeviations to minimize experimental error.\ncaption noted added: “Due to substrate inhibition induced by high concentrations of\nheme, the modified Michaelis-Menten equation incorporating an inhibition term was\nemployed for fitting”.\nIn the main text, we updated Fig. 3c and the description as follows:\n“During testing, decreased catalytic rates of HO and its variants were observed at\nhigher heme concentrations. Consequently, the Michaelis-Menten equation\nincorporating a substrate inhibition term was employed for fitting, yielding kinetic\nparameters including K.” (Line 173-177)\ni\nRegarding the original SI Fig. 6 mentioned in Reviewer's, severe substrate\ninhibition in RnBVR prevented observation of catalysis rate enhancement within the\ninstrumentally quantifiable substrate range. We have therefore supplemented the main\ntext with specific activity values based on the detectable maximum reaction rate (Fig.\n3c). Although these values are lower than the actual specific activity of RnBVR,\ncomparative analysis with HO and its mutants demonstrates that the BVR-mediated\nstep does not constitute the bottleneck in the overall cascaded process.\n4. Fig. 2A-C: these figures are not directly cited in the manuscript.\nResponse: Thanks a lot! The corresponding citation has been added at Line 137\n“We conducted phylogenetic analysis of HO (Fig. 2a), RP (Fig. 2b), and BVR (Fig.\n2c)”.\n5. Line 322: please delete the duplicated word (including)\nResponse: Thanks a lot! The duplicate word at Line 322 has been removed (now\nLine 334 in the revised manuscript).\n6. In the methods, could you provide how the lysate was made for the bioreactor?\nResponse: We thank the Reviewer for pointing this out! In our initial submission,\nwe included details for small-scale cell lysate preparation in the Supplementary\n7\ninformation but omitted the description for large-scale lysate preparation used during\ncatalysis.\nTo address this, we have added the following to the Methods section (Line 621-\n624)\n“Cell lysate for the 1 L scale was prepared using a high-pressure homogenizer\n(ATS AH-1500), operated at 850 bar with a flow rate of 5 L/h. The cell suspension was\npassed through the homogenizer three times to obtain the lysate, with chilled circulating\nwater maintained at 4 °C.”\nThe study sits at the intersection of synthetic biology, biocatalysis and enzyme\nengineering, where the authors reconstitute or redesign multi-enzyme pathways to\nmanufacture high-value natural products (the bile pigment bilirubin) from inexpensive\nprecursors. By combining five enzymes in a single, cell-free condition, the authors\nincrease bilirubin (a bile-acid precursor) titer from the long-standing ceiling of tens of\nmilligrams per liter to 1.68 g L⁻¹ in <3 h. This is a >20-fold jump that turns a\nniche diagnostic dye into a scalable valuable product. This breakthrough demonstrates how\npathway conditioning can unlock high-value tetrapyrrole manufacturing in synthetic\nbiology. The study couples physics-based modeling (DFT + MD) with experimentally-\ndriven cascade design. Quantum chemistry calculations detail how two coproducts of\nheme oxygenase—Fe²⁺ and CO—trigger biliverdin/bilirubin degradation or enzyme\nstalling; parallel bench work implements iron chelation, CO scavenging, and NADPH\nrecycling to avert those bottlenecks, delivering gram-scale output. In general, I consider\nthe study significant and innovative. Here are some comments:\n1. Relevant to the current study, a recent perspective focuses on the impact of physics-\nbased modeling (QM and MD) on enzyme engineering. The present work illustrates a\ncomplementary aspect—intermediate-guided engineering—showing that managing\nunstable or inhibitory pathway intermediates (Fe²⁺, CO) can be more decisive for yield\nthan enzyme optimization. I would suggest the author to conduct a super brief\ncomparison to the perspective by Jurich et al. (“Physics-based modeling in the new era\nof enzyme engineering,” Nat. Comput. Sci. 5, 279–291, 2025), highlighting this\ncontribution as a complementary design tier for synthetic biology.\nResponse: We sincerely thank Review for the suggestion to include a comparative\ndiscussion between this work and other enzyme engineering approaches. We agree that\nthis addition enhance the depth of the Discussion section.\nAccordingly, we have added the following statement to the revised manuscript\n(Line 503-511)\n9\n\"While the development of computational biology has established an effective\nparadigm for enhancing target protein properties (e.g., activity, stability) by employing\nphysical modeling to precisely describe atomic or electronic-scale dynamics and guide\nenzyme engineering55, this work demonstrates a complementary perspective in\nsynthetic biology. It highlights that in biocatalytic processes, intermediate-guided\nengineering—particularly in managing unstable or inhibitory pathway intermediates\n(such as Fe²⁺, CO, etc.)—can exert a more critical influence on yield improvement than\nenzyme optimization alone.\"\n2. Major conceptual critique: computation is confirmatory, not hypothesis-generating.\nBefore this study the community already knew that 1) free Fe²⁺ catalyzes oxidative\ncleavage of biliverdin/bilirubin, and chelators or acidification can suppress the damage;\n2) CO binds heme far more tightly than O₂, forming a photolabile heme-CO complex\nthat inhibits heme-oxygenase turnover; employing CO-scavenging enzymes is a\nrecognized workaround. It is true that the authors’ DFT and MD simulations further\nquantify these effects, but I would argue that these computations are primarily\nconfirmatory; the chelator (HEDP) and CODH choices appeared to be guided by prior\nempirical knowledge rather than by de-novo quantum chemistry insight. I would\nsuggest the authors to: 1) emphasize the design principle that the DFT results (homo-\nlumo gap and binding free energy) are uniquely poised to inform (some guidance by\nDFT that can not be obtained from simple empirical understanding); and/or 2) consider\nframing DFT more as a supportive approach to achieve the goal, rather than guiding the\ndevelopment of the approach.\nResponse: We sincerely thank the Reviewer's conceptual critique and their two\npivotal suggestions: emphasizing the design principles and positioning computational\nresults as supportive validation tools. We understand the concern that presenting\npreviously established chemical mechanisms as novel could reduce the significance of\nour work.\nRegarding Fe²⁺-mediated effects on bilirubin/biliverdin, studies remain limited.\nBilirubin can convert to biliverdin in the presence of Fe²⁺ and oxygen radicals\n(DOI: 10.1002/rcm.3822);\nExtensive research has addressed ROS-induced degradation of bilirubin/biliverdin\n(DOIs: 10.1016/j.freeradbiomed.2010.06.001; 10.1016/j.bmc.2015.11.012;\n10.1016/j.foodchem.2014.10.072), but these typically employ isolated systems without\nFe²⁺ involvement;\nSome reports note Fe²⁺-associated destabilization of bile pigments without\nmechanistic exploration (DOI: 10.1016/j.ab.2009.05.040)\nThus, our investigation of Fe²⁺-induced oxidative degradation provides some\nmechanistic insights to porphyrin chemistry, extending beyond confirmatory\nexperimentation. We respectfully request reconsideration of this work's novelty and\nsignificance.\nRegarding CO-heme binding, we acknowledge it is well-established. Previous\nstudies have documented evolutionary adaptations in heme oxygenase to mitigate CO\ninhibition, including steric constraints in the ligand-binding pocket and distal\nhydrophobic cavities for CO sequestration. Our DFT/MD simulations on CO inhibition\ntherefore serve primarily as confirmatory validation for the engineered rat heme\noxygenase mutant.\nWe implemented the Reviewer's suggestions through the following revisions:\nAbstract: Line 14, “two” → “the”; Line 15, deleted “and CO-inhibited heme\noxygenase activity”;\nMain text:\nLine 92-94, added “On the other hand, we provided additional theoretical insights\ninto the CO inhibition phenomenon, which is known in the community”;\nLine 94, deleted “We also demonstrated that”;\nLine 354, “CO-mediated inhibition of heme oxygenase as a hidden bottleneck in\ncascade biosynthesis” → “Effects of gaseous byproduct: verification of CO inhibition”;\nLine 361-363, “we hypothesized that CO accumulation may impede enzymatic\n11\nturnover” → “assessment of CO inhibitory effects on the established heme oxygenase\nsystem is required”;\nLine 377, “explore” → “verify”;\nLine 386-388, added “This geometric disparity induces steric hindrance consistent\nwith crystallographic observations, ultimately resulting in significantly lower affinity\nof HO for Heme-CO than for Heme-O₂.”, [DOI: 10.1021/bi027268i];\nLine 396, added prefix: “supplementary verification of”;\nLine 400, deleted: “thus establishing CO as a previously unrecognized bottleneck\nin cascade biosynthesis of BR”;\nLine 468, “two” → “the”;\nLine 470, deleted “and CO-induced inhibition of HO enzyme”;\nLine 488, “identification” → “verification”;\nLine 490-492, “its direct impact on HO catalysis had not been mechanistically\nexplored” → “we offer further insights into this inhibition mechanism through\nsimulation tools”.\nLine 495, deleted “novel and”;\nLine 521, “two” → “the”.\nThese revisions explicitly reframe the CO-related DFT/MD simulations as\nvalidatory tools supporting design guidance, rather than proposing a novel hypothesis.\nWe believe these amendments address your concerns and enhance the manuscript's\nclarity and quality.\n3. Other technical details. Spin-state coverage: only the quintet state is reported. It's\nimportant to consider the triplet state as it may lie near-degenerate for Fe²⁺ porphyrins.\nResponse: We sincerely thank the Reviewer's insightful comments regarding the\nelectronic structure of porphyrins. In our initial quantum chemical calculations, all\npossible spin multiplicities (singlet, triplet, and quintet) were explicitly for each\nstructure. The lowest energy state identified from this ordering was selected as the\nground state for further investigation.\nfollowing in the revised manuscript:\n1) Re-optimized structures and conducted vibrational analyses at the B3LYP-D3/6-\n311G(d,p) level.\n2) Supplemented the computational data with two additional levels of theory,\nresulting in three distinct computational levels:\nBP86-D3/def2-TZVP\nPBE0-D3/def2-TZVP\nDSD-PBEP86-D3/def2-TZVPP (double-hybrid functional)\nFollowing the Reviewer’s suggestions, we observed that, for certain porphyrin\nstructures, the triplet and quintet states were energetically very close under different\ncomputational conditions. This occasionally led to inconsistencies in the predicted\nground-state spin, validating the Reviewer's concern and addressing an oversight in our\ninitial submission.\nImportantly, despite variations in ground-state spin across computational level, all\nmethods consistently demonstrated that CO is a stronger ligand than O₂ for Heme,\nwhich remains the central focus of our study. The average computed binding free\nenergies across the three levels were: CO to Heme: -23.2 kcal/mol and O₂ to Heme: -\n6.9 kcal/mol.\nRelevant manuscript sections have been updated as follows:\nLine 366-374 (Paragraph Replacement):\nOriginal: “DFT calculations revealed that CO forms thermodynamically more\nstable complex with heme (ΔG = -40.8 kcal/mol) than O₂ (ΔG = -23.3 kcal/mol) (Fig.\n7d), suggesting that CO could outcompete O₂ for the coordination site.”\nRevised: “Employing three distinct levels of theory, including the double-hybrid\nfunctional DSD-PBEP86-D3, it was confirmed that CO forms a thermodynamically\nmore stable complex with Heme compared to O₂ (Fig. 7d). This supports CO's role as\na potent competitive agent capable of displacing O₂ from its coordination site. Although\nthe ground-state spin differed under various computational conditions, the average\n13\nbinding free energy of CO to Heme (-23.2 kcal/mol) was significantly stronger than\nthat of O₂ (-6.9 kcal/mol).”\nFigure 7d: Re-plotted to include the free energy data for all computational levels\nand all spin states, highlighted in yellow in the revised manuscript.\nLine 602-603, added the method description “For key species, calculations were\nadditionally performed using two DFT methods: PBE0/def2-TZVP and DSD-\nPBEP86/def2-TZVPP.”.\n4. Functional/basis set validation. Functional/basis set – B3LYP-D3(BJ)/6-31G(d,p) is\nused here. Suggest to consider other range-separated functional and also DLPNO-\nCCSD(T) single-points on key species. This would enhance reliability.\nResponse: We sincerely thank the Reviewer's for the insightful suggestions\nregarding computational chemistry. We recognized the relying on a single\nfunctional/basis set combination may introduce deviations from the actual situation. To\nstrengthen the credibility of our computational hypotheses, we attempted calculations\nat the DLPNO-CCSD(T)/def2-TZVPP level using the ORCA. However, due to the\nhardware limitations (48 CPU cores and 96 GB of RAM in total), a single complete\ncalculation would require over 7 days, exceeding our current capacity. As such, we were\nunable to obtain energy data for all key species within the revision timeframe. We plan\nto update our computational resources to perform these calculations in future work.\nTo address the Reviewer’s concerns in the meantime, we re-optimized the\ngeometries of all key species across all possible spin states at the B3LYP-D3/6-\n311G(d,p) level. Single-point energy calculations were then calculated using BP86-\nD3/def2-TZVP, PBE0-D3/def2-TZVP, and DSD-PBEP86-D3/def2-TZVPP. The\nupdated computational results have been incorporated into Figure 5, with the following\nmanuscript revisions:\nFig. 5c and Fig. 5d: Re-plotted with updated captions\nLine 252: “(-126.5 kcal/mol)” → “(average -132.7 kcal/mol, Fig. 5c)”.\nLine 253: “(-133.5 kcal/mol)” → “(average -153.4 kcal/mol, Fig. 5d)”.\n(deBV-Fe-O₂-OH) is -11.8 kcal/mol.\nLine 263 and Supplementary Fig. 15: Updated binding energy for species 15-5\n(deBR-Fe-O₂) is -26.1 kcal/mol.\nLine 263-270: Removed “(Fig. 5d)” and added:\n“To prevent potential biases arising from a single functional/basis set combination,\nenergy validation was performed across all spin states using three distinct\nfunctional/basis set combinations (Figs. 5c and 5d). While PBE0 calculations tend to\noverestimate Fe²⁺ coordination free energies and predict slightly different ground-state\nspins within these three DFT methods, all levels consistently indicate that coordination\nby deprotonated BV and BR, as well as subsequent axial ligand binding, is energetically\nfavorable and potentially spontaneous.”\nAlthough PBE0-D3 results are significantly higher than those from the other two\ncomputational levels, we have included them for transparency. We believe these\nrevisions address the Reviewer’s concerns regarding functional/basis set validation, and\nwe plan to complete DLPNO-CCSD(T) validation in future work.\nWe hope our updates satisfactorily address the Reviewer’s comments and look\nforward to further feedback.\n5. Limited sampling from MD. A single 200 ns replica of HO-1 is short for a membrane-\nassociated protein; multiple replicas (3 x 200 ns) or enhanced sampling would firm up\nthe Ser142 steric-block hypothesis.\nResponse: Thank Reviewer for highlighting the issue of insufficient sampling in\nour molecular dynamics (MD) simulations. To enhance the reliability of our results, we\nperformed three independent 200-ns MD simulations for both the Heme-O₂-HO and\nHeme-CO-HO complexes, maintaining consistency with our previous simulations in\nterms of force field, protonation states and solvent models.\nBased on the additional simulation results, we have updated Supplementary Fig.\n20 and its caption, as well as the description of MD method (Line 198-200) in the\nmmPBSA binding free energy calculations, resulting in 600 structural snapshots (200\n× 3). The binding free energy data were updated in Fig. 7e, and the corresponding\nmethod description revised (Line 202-210). The average binding free energy change\nafter triplicate simulations is now -2.50 kcal/mo for Heme-O₂, and 5.54 kcal/mol for\nHeme-CO, it is updated in Line 381. With the positive value for Heme-CO further\nsupporting its reduced binding affinity.\nResidue decomposition data from the binding free energy calculations were used\nfor protein coloring, updating in Fig. 7f and Fig. 7g. Fig. 7h now illustrates the\ncontributions of substrate-pocket residues to binding free energy, highlighting the\ninfluence of residues beyond S142 on ligand affinity. Corresponding textual\ndescriptions have been updated in Line 389-390 and in the caption of Fig. 7.\nFinally, the RMSD and Rg values from the three simulations were used to\nconstruct the free energy landscapes in Fig. 7i and Fig. 7j.\nIn summary, the addition of three independent 200-ns simulations provides\nstronger support for our conclusion regarding the reduced affinity of Heme-CO for the\nprotein.\n6. Data availability. In addition to the Cartesian coordinates, sharing input decks, force-\nfield parameters, and MD trajectories will enhance the reproducibility of the current\nwork.\nResponse: Thank Review very much! Providing essential input file information\nfor the QM calculations and molecular dynamics simulations will significantly enhance\nthe reproducibility and data availability of this work, which is critical for sharing high-\n sharing high-\nquality research within the scientific community.\nAccordingly, in the revised Supplementary information (pp. 45–111), we have\nincluded the net charge, spin multiplicity, and Cartesian coordinates of all molecules\nused in QM calculations, as referenced in both the main manuscript and Supplementary\ninformation. Each structure has been assigned a unique identifier (name or number) for\n16\nclarity. Detailed mappings can be found in Supplementary Figs. 14, 15, Fig. 17, and\nFig. 18, with captions refined for consistency.\nFor the molecular dynamics simulations, we will upload the force field parameter\nfiles for the ligands Heme-CO and Heme-O2 via the submission system; and given that\nthe trajectory files are substantial in size, they are available from the authors upon\nreasonable request.\nThese files will be made available for download on the Nature Communications\nwebsite upon publication, barring unforeseen issues. We believe these supplemental\ndata will satisfy the journal’s and reviewers’ requirements for data availability.\nReviewer #3:\nThe manuscript titled “Mechanistic Multi-Enzyme Engineering for High-Yield\nBilirubin Biosynthesis” submitted by Zhang et al., in the Nature Communications\nreports a comprehensive study achieving gram-scale bilirubin production through a\nrationally engineered multi-enzyme cascade. The identification of Fe²⁺-induced\noxidative degradation as a key bottleneck, supported by high-level computational\nanalysis and biochemical validation, represents a major conceptual advance. This\nmanuscript offers a systems-level approach to pathway optimization that could broadly\nbenefit synthetic biology and biocatalysis. However, several key areas require\nadditional clarification and validation before acceptance.\n1. The title emphasizes “Mechanistic Multi-Enzyme Engineering”, which\nsuggests extensive protein engineering (e.g., mutagenesis, directed evolution, or\nrational redesign of catalytic sites). However, the actual enzyme engineering performed\nin this work is primarily:\n• Truncation of hydrophobic regions of RnHO-1 and RnCPR to improve soluble\nexpression in E. coli, while retaining catalytic efficiency.\n• Cascade optimization, including tuning HO/CPR ratios and managing substrate\ninhibition in BVR.\n• System-level interventions, such as Fe²⁺ chelation, pH modulation, and cofactor\nrecycling to mitigate byproduct interference.\nWhile these are valuable and effective strategies, the novelty lies more in pathway\nand system-level engineering rather than in traditional protein engineering. To avoid\npotential misinterpretation, the authors might consider clarifying this distinction in the\ntitle, abstract, or introduction.\nResponse: We sincerely appreciate the reviewer for raising this important point\nregarding the accuracy of our original title. As the reviewer correctly notes, our work\nfocuses on the truncation mutagenesis of pathway enzymes, cascade construction and\noptimization, and system-level interventions such as managing byproduct interference.\nrevised the title from:\n“Mechanistic Multi-Enzyme Engineering for High-Yield Bilirubin Biosynthesis”\nto\n“Breaking Barriers in Bilirubin Biosynthesis via Systemic Engineering and\nGlobal Regulation”.\nWe believe the revised title more accurately conveys the emphasis of this study\nand avoids potential misinterpretation. We have also carefully reviewed the Abstract\nand Introduction and found no statements likely to mislead readers, but we remain\ngrateful for this constructive suggestion and will be glad to further adjust if the reviewer\nhave additional concerns.\n2. The manuscript aims to provide “additional theoretical insights into the CO\ninhibition phenomenon” and to “mechanistically explore its direct impact on HO\ncatalysis.” However, there is a crucial inconsistency between the mechanistic claims\nand the reported computational data:\n• The Discussion section states that “the stable coordination of CO to heme limits\nproper substrate positioning in the HO catalytic pocket, preventing O₂ binding and\nactivation.” This explanation assumes that the heme–CO complex binds favourably\nwithin the HO active site and thereby inhibits catalysis.\n• In contrast, the MD simulations report binding free energies of +5.54 kcal/mol\nfor the HO–heme–CO complex (unfavourable) and –2.50 kcal/mol for HO–heme–O₂\n(favourable). A positive binding free energy indicates that the CO-bound complex is\nunstable within the enzyme environment, directly contradicting the claim of “stable\ncoordination.”\n• In line 369 it is stated that, “it was confirmed that CO forms a\nthermodynamically more stable complex with Heme compared to O₂”\nThis internal inconsistency regarding a central mechanistic finding for a key\nbottleneck is a major weakness that fundamentally undermines the paper's scientific\n2\nconclusions.\nResponse: We sincerely thank the reviewer for this insightful comment. We\nrecognize that our original wording did not adequately convey of the relationship\nbetween the computational results and the mechanistic explanation, which may have\ncaused confusion.\nThe inhibitory impact of CO on the heme oxygenase (HO) catalysis operates\nthrough two main mechanisms:\n(1) In solution: CO binds more strongly than O₂ to free heme (Fig. 7d), forming a\nthermodynamically stable heme-CO complex. This shifts the equilibrium toward the\nCO-bound form.\n(2) Within HO: the pre-formed heme–CO complex is poorly accommodated in\nthe HO catalytic pocket, as reflected by the higher binding free energy compared to the\nheme–O complex (Fig. 7e). Even if binding occurs, the strong coordination of CO at\n2\nthe axial iron position prevents O displacement and subsequent oxidative activation.\nbinding in HO, leading to catalytic inhibition.\nTo eliminate this inconsistency, we have revised the Discussion section as follows\n(Line 405-413):\nOriginal text: “the stable coordination of CO to heme limits proper substrate\npositioning in the HO catalytic pocket, preventing O₂ binding and activation. This\ncompetitive inhibition is further amplified by CO accumulation in solution, where it\nacts as a heme-sequestering agent.”\nRevised text: “CO acts as a stronger ligand than O₂, forming a stable heme–CO\ncomplex with free heme in the solution. However, pre-formed complex exhibits steric\nincompatibility with the HO catalytic pocket. Even if accommodated, CO occupies the\naxial position of the iron center, where its strong coordination prevents O₂ displacement\nand activation, thereby halting the subsequent catalytic cycle.”\nWe believe this revised description resolves the apparent inconsistency and more\naccurately reflects the computational and mechanistic data.\n3\n3. Ground-State Spin States: The authors mention that \"PBE0\ncalculations...predict slightly different ground-state spins within these three DFT\nmethods\" for Fe²⁺-containing systems. However, they do not elaborate on which spin\nstate was ultimately chosen for the reported energies, how these differences were\nreconciled, or what implications these variations have for the reliability of their\nelectronic structure calculations. For metalloenzyme systems, spin state can\nsignificantly influence reactivity, and this lack of detailed explanation weakens\nconfidence in the computational methodology.\nResponse: We sincerely thank the Reviewer for their expert comments regarding\nthe spin-state calculations in Fe2+ systems. We agree that spin states are critical factor\nin metalloenzyme reactivity and that accurate prediction remains a longstanding\nchallenge in computational chemistry. Benchmark studies (e.g., Radon et al.\nhttps://doi.org/10.1039/d4cp01327a) have shown that DFT functionals often introduce\nsystematic biases, with some over-stabilizing intermediate-spin states (e.g., B3LYP-D3,\nTPSSh-D3) and others favoring high-spin states (e.g., M06L-D3).\nGiven the lack of experimental reference data and the impracticality of higher-\nlevel methods such as CCSD(T), we employed a multi-method strategy. Specifically,\nthree independent DFT functional/basis set combinations were used to compute\nground-state energies and corresponding binding free energies. The final values\nreported represent averages across these methods, thereby reducing the likelihood of\nsystematic bias.\nImportantly, although the predicted ground states varied among methods (e.g.,\ndeBV–Fe was consistently quintet across all, whereas deBR–Fe showed quintet with\nBP86-D3 but triplet with PBE0-D3 and DSD-PBEP86-D3), all three approaches\nyielded the same qualitative outcome: Fe²⁺ coordination by deprotonated biliverdin and\nbilirubin, as well as subsequent O₂ binding, is energetically favorable. This consistency\nacross different ground-state descriptions supports the mechanistic plausibility of the\nproposed role of Fe²⁺.\n Lines 266–276: Added explicit description of the spin-state outcomes\nobtained with each functional.\n Lines 276–284: Reworded to emphasize that, despite spin-state variations,\nall methods converge on the same mechanistic conclusion.\n Lines 268 and 766–768: Added a reference to benchmark studies on spin-state\nprediction in metalloporphyrins.\nWe believe these revisions strengthen both the transparency and reliability of our\ncomputational analysis.\n4. Fundamental Flaw in Enzyme Kinetic Characterization (RnBVR): The article\nexplicitly states a critical limitation in characterizing one of the core enzymes in the\nbilirubin biosynthesis cascade. It mentions, \"strong substrate inhibition in the\ndownstream reaction catalyzed by RnBVR prevented accurate kinetic modeling\n(Supplementary Fig. 7), underscoring the importance of enzyme-specific optimization\nstrategies in cascade system design\". Furthermore, it clarifies, \"Consequently, kinetic\nparameters of RnBVR could not be reliably determined, and enzymatic activity was\nonly calculated at the lowest detectable substrate concentration permitted by the\ninstrument\".\nResponse: We sincerely thank the Reviewer for highlighting the importance of\naccurately characterizing RnBVR kinetics. We fully agree that this enzyme, which\ncatalyzes the biliverdin-to- bilirubin conversion, plays a central role in the cascade.\nAs mentioned in the manuscript, our current instrumentation limits biliverdin\nquantification to concentrations ≥50 μM. Below this threshold, noise from UV\ndetection precludes reliable measurement. Unfortunately, within the measurable\nconcentration range, substrate inhibition obscures the Michaelis–Menten relationship,\npreventing accurate determination of kinetic parameters.\nTo provide at least a comparative reference, we reported specific activity at the\nlowest detectable concentration (50 μM). Although this underestimates, the RnBVR\n5\nmaintains substantially higher activity than RnHO-1, supporting its effective function\nwithin the cascade.\nWe acknowledge that this dataset is not perfect and are actively developing more\nsensitive detection strategies to enable full kinetic characterization in future work.\nHowever, we believe the current results still provide a meaningful and sufficient robust\nbasis for the cascade design and conclusion presented here.\n5. Abstract (line 20) - “block O2-activation for accessing heme oxygenase” →\nsuggested “blocking O₂ activation at heme oxygenase.”\nResponse: Thank the comment a lot! We revised “block O -activation for\n2\naccessing heme oxygenase” to “blocking O₂ activation at heme oxygenase” (Line\n22).\n6. Discussion (line 501) - “Comparation” → should be “Comparison.”\nResponse: Thank this comment a lot! We revised “ Comparation ” to\n“Comparison” (Line 512).\n7. “the role of byproduct-metal interactions” (line 478) could be rephrased for\nclarity.\nResponse: We thank the reviewer for this helpful suggestion! The phrase “We\nuncovered that the role of byproduct-metal interactions” has been revised to “We have\nrevealed that the fundamental reason lies in the destabilization and degradation caused\nby the interaction between the product and the metal byproduct” (Line 488-491).\n8. Line 12: coverts → converts\nResponse: We revised “coverts” to “converts” (Line 14).\n9. Line 13: Systematically enzyme should likely be Systematically, enzyme\nResponse: Thanks a lot! We revised “Systematically enzyme” to “Systematically,\nenzyme” (Line 15).\n6\n10. Line 524: nature products → natural products\nResponse: We revised “nature products” to “natural products” (Line 535).\n11. Line 176 in SI: pymol → PyMol\nResponse: We revised “pymol” to “PyMol” (Line 179) in SI.\n12. Line 179 in SI: Molecular dynamic simulation should be Molecular dynamics\nsimulation\nResponse: We revised “Molecular dynamic simulation” to “Molecular dynamics\nsimulation” (Line 182) in SI.\n13. Line 436: fully integrate should be fully integrated\nResponse: We revised “fully integrate” to “fully integrated” (Line 448).\n14. Line 816: conducted DFT calculation should be conducted DFT calculations\nor conducted a DFT calculation\nResponse: Thank Review very much! We revised “conducted DFT calculation”\nto “conducted DFT calculations” (Line 830).\n7"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed my concerns.", "reviewer_2": "The authors have adequately addressed my review comments from last round of review. It's understandable that energetic calculations with DLPNO-CCSD(T) would take substantial amount of time, and alternative strategies employed by the authors (examination with multiple other DFT methods) are acceptable. I no longer have any further comments on this matter.", "reviewer_3": "The manuscript titled “Mechanistic Multi-Enzyme Engineering for High-Yield Bilirubin Biosynthesis” submitted by Zhang et al., in the Nature Communications reports a comprehensive study achieving gram-scale bilirubin production through a rationally engineered multi-enzyme cascade. The identification of Fe²⁺-induced oxidative degradation as a key bottleneck, supported by high-level computational analysis and biochemical validation, represents a major conceptual advance. This manuscript offers a systems-level approach to pathway optimization that could broadly benefit synthetic biology and biocatalysis. However, several key areas require additional clarification and validation before acceptance.\n1. The title emphasizes “Mechanistic Multi-Enzyme Engineering”, which suggests extensive protein engineering (e.g., mutagenesis, directed evolution, or rational redesign of catalytic sites). However, the actual enzyme engineering performed in this work is primarily:\n• Truncation of hydrophobic regions of RnHO-1 and RnCPR to improve soluble expression in E. coli, while retaining catalytic efficiency.\n• Cascade optimization, including tuning HO/CPR ratios and managing substrate inhibition in BVR.\n• System-level interventions, such as Fe²⁺ chelation, pH modulation, and cofactor recycling to mitigate byproduct interference.\nWhile these are valuable and effective strategies, the novelty lies more in pathway and system-level engineering rather than in traditional protein engineering. To avoid potential misinterpretation, the authors might consider clarifying this distinction in the title, abstract, or introduction.\n2. The manuscript aims to provide “additional theoretical insights into the CO inhibition phenomenon” and to “mechanistically explore its direct impact on HO catalysis.” However, there is a crucial inconsistency between the mechanistic claims and the reported computational data:\n• The Discussion section states that “the stable coordination of CO to heme limits proper substrate positioning in the HO catalytic pocket, preventing O₂ binding and activation.” This explanation assumes that the heme–CO complex binds favourably within the HO active site and thereby inhibits catalysis.\n• In contrast, the MD simulations report binding free energies of +5.54 kcal/mol for the HO–heme–CO complex (unfavourable) and –2.50 kcal/mol for HO–heme–O₂ (favourable). A positive binding free energy indicates that the CO-bound complex is unstable within the enzyme environment, directly contradicting the claim of “stable coordination.”\n• In line 369 it is stated that, “it was confirmed that CO forms a thermodynamically more stable complex with Heme compared to O₂”\nThis internal inconsistency regarding a central mechanistic finding for a key bottleneck is a major weakness that fundamentally undermines the paper's scientific conclusions.\n3. Ground-State Spin States: The authors mention that \"PBE0 calculations...predict slightly different ground-state spins within these three DFT methods\" for Fe²⁺-containing systems. However, they do not elaborate on which spin state was ultimately chosen for the reported energies, how these differences were reconciled, or what implications these variations have for the reliability of their electronic structure calculations. For metalloenzyme systems, spin state can significantly influence reactivity, and this lack of detailed explanation weakens confidence in the computational methodology.\n4. Fundamental Flaw in Enzyme Kinetic Characterization (RnBVR): The article explicitly states a critical limitation in characterizing one of the core enzymes in the bilirubin biosynthesis cascade. It mentions, \"strong substrate inhibition in the downstream reaction catalyzed by RnBVR prevented accurate kinetic modeling (Supplementary Fig. 7), underscoring the importance of enzyme-specific optimization strategies in cascade system design\". Furthermore, it clarifies, \"Consequently, kinetic parameters of RnBVR could not be reliably determined, and enzymatic activity was only calculated at the lowest detectable substrate concentration permitted by the instrument\".\n5. Abstract (line 20) - “block O2-activation for accessing heme oxygenase” → suggested “blocking O₂ activation at heme oxygenase.”\n6. Discussion (line 501) - “Comparation” → should be “Comparison.”\n7. “the role of byproduct-metal interactions” (line 478) could be rephrased for clarity.\n8. Line 12: coverts → converts\n9. Line 13: Systematically enzyme should likely be Systematically, enzyme\n10. Line 524: nature products → natural products\n11. Line 176 in SI: pymol → PyMol\n12. Line 179 in SI: Molecular dynamic simulation should be Molecular dynamics simulation\n13. Line 436: fully integrate should be fully integrated\n14. Line 816: conducted DFT calculation should be conducted DFT calculations or conducted a DFT calculation"}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_3": "I have reviewed the revised manuscript and the authors have addressed all the comments. I have no further suggestions."}, "author_response": ""}]

Decision letter