Engineering analysis of multienzyme cascade reactions for 3ʹ‐sialyllactose synthesis

Abstract Sialo‐oligosaccharides are important products of emerging biotechnology for complex carbohydrates as nutritional ingredients. Cascade bio‐catalysis is central to the development of sialo‐oligosaccharide production systems, based on isolated enzymes or whole cells. Multienzyme transformations have been established for sialo‐oligosaccharide synthesis from expedient substrates, but systematic engineering analysis for the optimization of such transformations is lacking. Here, we show a mathematical modeling‐guided approach to 3ʹ‐sialyllactose (3SL) synthesis from N‐acetyl‐ d‐neuraminic acid (Neu5Ac) and lactose in the presence of cytidine 5ʹ‐triphosphate, via the reactions of cytidine 5ʹ‐monophosphate‐Neu5Ac synthetase and α2,3‐sialyltransferase. The Neu5Ac was synthesized in situ from N‐acetyl‐ d‐mannosamine using the reversible reaction with pyruvate by Neu5Ac lyase or the effectively irreversible reaction with phosphoenolpyruvate by Neu5Ac synthase. We show through comprehensive time‐course study by experiment and modeling that, due to kinetic rather than thermodynamic advantages of the synthase reaction, the 3SL yield was increased (up to 75%; 10.4 g/L) and the initial productivity doubled (15 g/L/h), compared with synthesis based on the lyase reaction. We further show model‐based optimization to minimize the total loading of protein (saving: up to 43%) while maintaining a suitable ratio of the individual enzyme activities to achieve 3SL target yield (61%–75%; 7–10 g/L) and overall productivity (3–5 g/L/h). Collectively, our results reveal the principal factors of enzyme cascade efficiency for 3SL synthesis and highlight the important role of engineering analysis to make multienzyme‐catalyzed transformations fit for oligosaccharide production.


| INTRODUCTION
Sialo-oligosaccharides have received increased attention as healthpromoting ingredients for food and feed use (Faijes et al., 2019;Lu et al., 2021). Their synthesis as industrial products (e.g., 3ʹ-sialyllactose, 3SL) drives an emerging biotechnology for the mass production of structurally defined, complex carbohydrates (Bode et al., 2016). 3SL is one of simplest of the human milk oligosaccharides (X.  and is currently considered strongly for commercial use in infant formula (Bych et al., 2019). Core task of every sialo-oligosaccharide process is to provide from expedient substrates the naturally scarce sialic acid (e.g., N-acetyl-D-neuraminic acid, Neu5Ac) in a form usable for enzymatic glycosylation (Fessner, 2015). The glycosylation involves stereo-and regioselective transfer of the sialic acid to the acceptor oligosaccharide (e.g., lactose) (R. Na et al., 2021;Weijers et al., 2008).
It is typically catalyzed by a sialyltransferase . The enzyme uses a cytidine 5ʹ-monophosphate (CMP)-activated sialic acid (e.g., CMP-Neu5Ac) as the substrate which for reason of process cost effectiveness cannot be added as a reagent, but must be synthesized directly in the reaction . Bio-catalysis in multienzyme cascades, integrating glycosylation with supply of CMP-Neu5Ac, is therefore central to the development of sialo-oligosaccharide production systems, irrespective of whether isolated enzymes or live whole cells are used (Schelch et al., 2020). Multienzyme transformations have been established for sialo-oligosaccharide synthesis with both systems (Schelch et al., 2020), but systematic engineering analysis of the reaction cascades for the optimization of such transformations is lacking. Kinetic modeling-based approaches are important engineering tools for coping with the inherent complexity of enzymatic cascades for efficient process development (Kitamura et al., 2020;Schmideder et al., 2017;Xue & Woodley, 2012;Zhong et al., 2017;Zimmermann et al., 2007). However, use of kinetic modeling for studying the assembly of oligosaccharides by sugar nucleotide-dependent transferases is generally scarce (Mahour et al., 2018;Rexer et al., 2017) and is missing entirely with sialo-oligosaccharides. Fundamental engineering problems, like the identification of main bottlenecks of conversion efficiency of sialooligosaccharide-producing enzyme cascade reactions, thus remain largely unaddressed. Relevant progress is important to advance the biocatalytic synthesis by using these enzymatic systems. It furthermore provides essential guidance to metabolic engineering efforts in cell factory development (e.g., 3SL-producing Escherichia coli; Faijes et al., 2019;Lu et al., 2021).
While the presence of acceptor substrate reduces this hydrolase activity (Schmölzer et al., 2017), it remains a significant side reaction. PdST can also act as a sialidase. The sialidase reaction, which overall releases Neu5Ac and lactose from 3SL, requires CMP (Mehr & Withers, 2016). It is believed to proceed via intermediary CMP-Neu5Ac formed through the reverse transferase reaction (3SL + CMP → CMP-Neu5Ac + lactose) (Mehr & Withers, 2016). The CMP-Neu5Ac is then hydrolyzed. The sialidase activity of PdST is found low when assayed directly from 3SL (Schmölzer et al., 2013) but its possible role in a synthesis reaction cannot be excluded.
In this study, we performed systematic engineering analysis for 3SL synthesis by NAL/CSS/PdST and SiaC/CSS/PdST cascade reactions. Through quantification of all relevant reactants from the enzymatic conversions, we obtained detailed time-course data for the multistep transformations and used them for kinetic model development. A mass action-based, Michaelis-Menten type kinetic model (Bulik et al., 2009) was used to describe the complete "reaction dynamics" in 3SL syntheses by the two enzymatic systems. Due to kinetic advantages (K M for ManNAc) of the SiaC compared with the NAL reaction, the 3SL productivity was generally enhanced when the SiaC/CSS/PdST system was used. The validated Michaelis-Menten model was employed for facile reaction "optimization" through an innovative in silico approach. A constrained random simulation of a large set of reaction conditions (≥10 5 ) with variable enzyme ratios was carried out. From the results of this high-throughput computational sampling, as a typical objective for the optimization of enzyme cascade reactions, conditions were selected that involved the least usage of total enzyme to satisfy a given conversion task. The modelpredicted optima were verified experimentally. Collectively, our results demonstrate kinetic modeling-guided development of efficient enzyme cascades for 3SL synthesis. Our findings highlight the important role of such analysis in the optimization of multienzyme biocatalysts for oligosaccharide production.

| Expression and purification of enzymes
PdST was obtained as described by Schmölzer et al. (2013). Expression was done similarly for all enzymes (1 mM isopropyl β-D-thiogalactopyranoside; 20 h) except that 25°C (SiaC, NAL, and ManNAcDH) or 37°C (CSS) was used during induction. Purification was done by His-tag affinity chromatography (Online Supporting Information). Enzyme purity was verified by sodium dodecyl sulfate-polyacrylamide gel electrophoresis ( Figure S1). Protein was determined with Roti-Quant reagent (Roth) referenced to bovine serum albumin (BSA). Enzyme stock solutions were stored (~20 mg/ml; buffer) at −20°C without loss of activity for at least 21 days.

| Activity assays
Assays were conducted in duplicate in 200 µl total volume of 100 mM Tris/HCl buffer, pH 8.0, containing 1 mg/ml BSA. Temperature (37°C) and agitation rate (450 rpm) were controlled in a Thermomixer comfort

| Enzymatic assay for ManNAc
This was adapted from Klermund et al. (2016) and performed in 96well microtiter plates. Sample (20 µl) from NAL or SiaC reaction was added to 145 µl of 100 mM Tris/HCl buffer (pH 8.0). NAD + (1 mM; 20 µl) and ManNAcDH (20 µl; 0.7 mg/ml) were added and incubation done for 20 min at 30°C and 450 rpm. NADH was measured at 340 nm using a multimode microplate reader (BMG Labtech FLUOstar Omega). The ManNAcDH forms one NADH for each ManNAc oxidized. Activity of NAL or SiaC was determined from the ManNAc consumed over time.
In Michaelis-Menten models, the two K M parameters associated with two substrate reactions were assumed to be independent one from another.  (Schmölzer et al., 2013). At low lactose concentrations at around the K M (≤1.5 mM), the ratio between hydrolysis and transfer (R h ) was approximately 0.17 (Schmölzer et al., 2015). To model the PdST reaction, we therefore added R h and K eq as additional fit parameters. The R h was assumed constant (not variable with concentration of acceptor substrate). The low "sialidase" activity of PdST ( conditions between assays and synthesis can affect all parameters (V max , K eq , and K M ), it is mostly the V max that governs the conversion rate in a cascade (see the Section 3.3 on the K eq of NAL later).
The fitting used MATLAB's lsqnonlin (nonlinear least-square method with trust-region-reflective algorithm) utility. In total, 5000 independent fits were performed. Each fit started by selecting random starting values for each parameter from individually predefined sets, containing 50 values and designed as follows. V max starting values were evenly distributed in the range ±50% of the V max from literature or the V max determined with the herein used assays (Table 1). R h starting values were evenly distributed between 0 and 0.1 (Schmölzer et al., 2015). Starting values for K eq were logarithmically distributed over three orders of magnitude for the NAL (range 10 −1 to 10 1 ) and the PdST reaction (range 10 0 -10 2 ).
Iterative fitting continued until the step size of the objective function fell below step size tolerance (MATLAB default: 10 −6 ). The number of used iterations never reached the MATLAB default maximum of 400.
The top 5% of the obtained fit parameter combinations based on residual error were further analyzed. After check for plausible agreement with the experimentally determined parameter (Table 1), the median of each series of estimates was used for optimization.

| Model guided optimization
The objective was to minimize the total protein loading from a sui-   Table S2.

| Enzyme characterization for use in cascade reactions
Essential requirement for the planned engineering analysis was that the enzyme preparations used were well defined and characterized.

SCHELCH ET AL.
| 4295 CMP-Neu5Ac (the most labile among the reactants of the cascade reaction) was best in the pH range 8.0-11.00 (Beau et al., 1984), we chose a "working pH" of 8.0. The CSS reaction requires Mg 2+ . We added 20 mM MgCl 2 in combination with 25 mM CTP as the substrate. A higher concentration of Mg 2+ (tested up to 40 mM with 25 mM CTP) neither enhanced the enzymatic rate nor did it improve the CMP-Neu5Ac yield. We also showed that in coupled reactions of SiaC (20 mM PEP) and CSS (25 mM CTP) no more than 20 mM Mg 2+ were necessary for maximum conversion rate and product yield.
L-Cysteine (0.2 mM) was added for CSS activity and stability. In our hands, the L-cysteine was equally effective as the dithiothreitol used in the literature (Mizanur & Pohl, 2008 As reference for optimization with each cascade system, a "common sense" synthesis reaction was designed in which equal volumetric activities (0.6 U/ml) of all enzymes were used. ManNAc and lactose were supplied at 20 mM each. PYR was used at 50 mM, following notion from earlier studies that 2.5-fold molar excess of PYR over ManNAc can drive the Neu5Ac formation . CTP was used in slight excess over ManNAc (25 mM), to account for CMP-Neu5Ac hydrolyzed by the PdST (Schmölzer et al., 2013). We estimated that the overall conversion (ManNAc → 3SL) proceeding at maximum rate (0.6 mM/min) in each enzymatic step would require approximately 100 min (= 20 × 3/0.6) to complete. Reactions were analyzed for 4 h.

| Reactant quantification for time-course analysis
Modeling of cascade reactions is best supported by analytical quantification of all involved, process-relevant reactants. Despite the widespread use of the NAL/CSS/PdST cascade reaction for sialooligosaccharide synthesis (Malekan et al., 2013;Tasnima et al., 2019;Yu et al., 2009Yu et al., , 2011, the reaction analysis was typically restricted to product, acceptor substrate and sometimes both (Lau et al., 2011;. Intermediates were not determined. We here therefore set out to quantify ManNAc, Neu5Ac, CTP/CMP, CMP-Neu5Ac, and 3SL. Based on close mass balance confirmed experimentally, PYR, PEP, lactose, and pyrophosphate were considered redundant and not measured routinely. Using TLC for preliminary assessment of the conversion ( Figure S2), we developed HPLC analytical protocols to determine CMP-Neu5Ac and nucleotides separate from the carbohydrates. In each protocol, a dedicated sample preparation was important to optimally suit the subsequent HPLC analysis. Ion-pairing on reversed C-18 stationary phase (nucleotides) or ligand exchange (carbohydrates) was used for separation ( Figure 2). To avoid overlap between CTP and 3SL in elution ( Figure S3), the nucleotides were hydrolyzed by a phosphatase before carbohydrate analysis. We noted that due to baseline shift caused, the acetonitrile used for enzyme inactivation interfered with the carbohydrate analysis. Heat inactivation was used instead. The analytical procedures were applicable to NAL/CSS/PdST as well as SiaC/CSS/PdST cascade reactions.  Table S4.

| SiaC and NAL cascade reactions for 3SL synthesis
The NAL reaction (Figure 3c,d) was overall slower (∼2-fold) than the SiaC reaction. The 3SL release rates in the initial 20 min were 7 g/(L h) and 15 g/(L h), respectively. The 3SL time-course featured a substantial decrease in the production rate already at low degrees of substrate conversion (≤ 50%; 45-60 min). Neu5Ac did not accumulate and CMP-Neu5Ac did so only slightly. From these reactant profiles, the NAL reaction appeared to have been limiting overall.
To exclude that enzyme inactivation could have restricted the attainable degree of substrate conversion and the 3SL yield, we measured the stability of each enzyme under the bulk conditions of the reaction in the absence of substrates. As shown in Table S3, activity loss was negligible over 1 h. Additionally, we analyzed samples taken directly from the reaction after 2 h. There was only a minor decrease in the activity (≤ 20%) of the enzymes (Table S5).

| Modeling of the enzymatic cascade reactions
We used mass-action kinetics with parameters (K M ) accounting for dependence of the reaction rate on the substrate concentration (Equation 1). As shown in Figure 3, using K M values from literature (Table 1) (Table 1). The K eq for the transferase reaction of PdST ( Figure 4e) was in a defined narrow range far on the product side.
Its average value of 158 was in accordance with literature (Schmölzer et al., 2013), reporting a ratio of approximately 300 for the enzymatic rates of forward (24 s −1 ) and reverse sialyl transfer (=sialidase activity in the presence of CMP; 0.08 s −1 ) from CMP-Neu5Ac to lactose. The estimated hydrolysis/transfer coefficient R h ( Figure 4f) was consistent with earlier results of initial rate analysis (Schmölzer et al., 2015). In contrast to the other parameters ( Figure 4), the K eq for the NAL reaction (Figure 4g) was not well defined as estimated from fitting. This apparent issue was resolved by understanding that in a cascade reaction, the K eq on an intermediate step can be neglected if a constant removal of the product is ensured (Ricca et al., 2011). Figure S4 shows simulated time courses of ManNAc consumption by a hypothetical NAL reaction in which K eq and K M for ManNAc were variable. The results reveal that variation in the K eq was without influence, explainable on account of the "pull" from the effectively irreversible CSS reaction. By contrast, lowering the K M for ManNAc resulted in a significant increase the ManNAc consumption rate ( Figure S4). The fitting results additionally revealed that in both cascade transformations but especially when SiaC was used (Figure 3a,b), the overall sialidase reaction, that is, the reverse PdST reaction coupled to hydrolysis of the CMP-Neu5Ac thus formed, contributed to limitation of the 3SL yield after 2 h.

| Modeling-based optimization of the enzyme loading
Objective for the model-based optimization was minimized total amount of protein used to synthesize 3SL as in the reference experiment (±10% tolerance) in a 2-h reaction. Thus, the TTN (g 3SL/g total protein) would be maximized. Concretely, the 3SL  Table S2.
F I G U R E 4 Boxplots of the fitted parameters derived from the top 5% simulations based on the sum of residual errors. The median is indicated by a black line, while the mean is shown in color and boxes extend from the 25th to the 75th percentile of each group's distribution. Whiskers show the 10th and the 90th percentile, respectively. The 5th and 95th percentiles are plotted as red dots. CSS, cytidine 5ʹmonophosphate-sialic acid synthetase; NAL, Neu5Ac lyase; Neu5Ac, N-acetyl-D-neuraminic acid; PdST, a2,3-sialyltransferase from P. dagmatis; SiaC, sialic acid synthase   simulations. The exact conditions are summarized in Table S2 and the results are shown in Figure S5. The excellent agreement between measured and simulated time-course data for both cascade reactions was strong evidence for model verification.

| CONCLUSIONS
Modeling-based approach to the optimization of three-enzyme cascade reactions for 3SL synthesis was presented. Besides the NAL cascade known from earlier studies (Schelch et al., 2020), the SiaC cascade was used here for the first time in bio-transformations in vitro. The optimization strategy was innovative: it was built on timecourse simulations with a mass action-controlled kinetic model that were done in substantial computational bulk (≥10 5 conditions) to screen a large, rationally defined operational space. The simulation results obtained thus enabled important optimization tasks to be analyzed flexibly. This was demonstrated for maximized TTN, based on individually optimized enzyme loadings, to reach a predefined conversion target. The optimized SiaC cascade gave higher 3SL yields (79% compared with 65%) and productivity (2-fold) than the optimized NAL cascade, and its corresponding TTN was almost three-fold higher. Comparison can be done with literature at the level of similar product concentration formed (Table S5). This shows that the computationally identified and experimentally verified conditions of optimized 3SL synthesis represented improvement by at least one magnitude order for the TTN, the productivity or both. The engineering analysis shown here can be generally relevant to promote the field of systems bio-catalysis (Fessner, 2015;France et al., 2017;Schmidt-Dannert & Lopez-Gallego, 2016;Yang et al., 2019), working with multienzyme cascade reactions in vitro (Lau et al., 2011;Malekan et al., 2013;Tasnima et al., 2019; but also in context of whole-cell metabolism (Faijes et al., 2019;Lu et al., 2021;Sprenger et al., 2017). It can be important to unlock the full potential of glycosyltransferase cascade reactions Mestrom et al., 2019;Nidetzky et al., 2018;Schelch et al., 2020; for efficient use in oligosaccharide and glycoside production. Lastly, it can support the making of fundamental choices in the development of enzyme cascade transformations, in particular which reactions should be telescoped in one pot and which rather not (e.g., Klermund et al., 2017;Rexer et al., 2020). The NAL cascade provides an interesting example for it could involve synthesis of Neu5Ac (Kragl et al., 1991;Lin et al., 2013;Lv et al., 2017;Schmideder et al., 2017;Tao et al., 2011) spatiotemporally separated from the sialoside formation, or integrated with it as shown here.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.