Escherichia coli strain Top10 was used for routine plasmid DNA preparation. Aspergillus niger strains are all derived from ATCC11414 from the American Type Culture Collection (Rockville, MD, United States) and were grown on complete medium (CM) or potato dextrose agar (PDA) plates (JW and LL, 1991) at 30 °C for culture maintenance and spore preparation. About 5 × 10⁴ to 5 × 10⁵ spores were inoculated on CM agar (Petri dish) plates and incubated for 4 days at 30 °C. Spores were harvested by washing with 5–10 mL sterile 0.4% Tween H₂O and pelleted by centrifugation at 2,500 g for 5 min. The spores were re-suspended in the sterile 0.4% Tween H₂O and enumerated with a hemocytometer. The spore suspensions were used for agar-plate or liquid cultures. The shake flask cultures were performed at 30 °C, 200 RPM in a New Brunswick Innova 44R stackable incubator shaker (Eppendorf, Endfield, CT, United States) with Pyrex 125 mL or 250 mL glass Erlenmeyer flasks which were prepared by filling with 5% Contrad 70 (Decon Labs, Inc., King of Prussia, PA, United States) and soaked overnight to remove any potential residues on the inside surface of flasks prior to general dishwashing. Silicon sponge closures were used for all flask cultures. For 3-HP production, sterile modified production medium B (mRDM; 100 g/L glucose, 5 g/L (NH₄)₂SO₄, 0.11 g/L KH₂PO₄, 2.08 g/L MgSO₄.7H₂O, 0.13 g/L CaCl₂.2H₂O, 0.074 g/L NaCl, 4 mg/L CuSO₄.5H₂O, 110 mg/L FeSO₄.7H₂O, 14 mg/L MnCl₂.4H₂O, 26 mg/L ZnSO₄.7H₂O) (Dai et al., 2023) was used for cultivations.

Strains were constructed as described previously (Dai et al., 2023; Pomraning et al., 2021) with additional strains containing intermediatory or combinatorial expression of transgene vectors included here to provide diversity in expression of 3-HP pathway components. In this study, all transgene expression cassettes were prepared with Gibson assembly master mix (NEB, Ipswich, MA, United States) and the DNA fragments were isolated by PCR with Phusion high-fidelity DNA polymerase (Thermo Fisher Scientific, Waltham, MA, United States). Strain ABF_008348, the highest performing strain in the reference dataset was used as the host for validation of target genes. In this background, targets for overexpression (purU, mdhA, pyc2, and aat1) were constructed using the several strong constitutive promoters (cox1, mbfA, tef1, and ubi4) and integrated into the genome. Strains containing deletions of adhD, suclg1, and suclg2 were constructed using CRISPR-Cas9 system for single gene editing in A. niger (Yuan et al., 2024) to create large indels in the genes targeted for deletion. Protoplast preparation and chemical-mediated transformation followed the method described previously (Dai et al., 2013). All strains used in this study are shown in Table 1. Plasmid and strain construction is described in further detail in Additional File 1.

For quantification of extracellular metabolites in the spent medium samples were analyzed using high performance liquid chromatography (HPLC) equipped with a Waters 2,414 refractive index detector. A Bio-Rad Aminex HPX-87H ion exclusion column (300 mm × 7.8 mm), heated to 65C was used for analyte separation. Sulfuric acid (0.0045 M) was used as eluent at a flow rate of 0.55 mL/min. Intracellular metabolite extracts were completely dried under vacuum to remove moisture and chemically derivatized. Briefly, the extracted metabolites were derivatized by methoxyamination and trimethylsilyation (TMS), then the samples were analyzed by GC-MS. Sample preparation, instrument acquisition and data analysis was performed as previously reported (Pomraning et al., 2021; Kim et al., 2015). Global and targeted proteomics was performed as in previously established LC-MS/MS methods and data analysis workflows (Pomraning et al., 2021), except slight adjustments to the mass spectrometry acquisition settings in global proteomics. In global proteomics, peptide digests were analyzed using a Q Exactive Plus mass spectrometer (Thermo Fisher Scientific) in data-dependent acquisition mode. Mass spectrometer settings were as following: full MS (AGC, 3 × 10⁶; resolution, 70,000; m/z range, 300–1800; maximum ion time, 20 m); MS/MS (AGC, 1 × 10⁵; resolution, 17,500; m/z range, 200–2000; maximum ion time, 50 m; TopN, 12; isolation width, 1.5 Da; dynamic exclusion, 30.0 s; collision energy, NCE 30).

BMCA was performed following established protocols (McNaughton et al., 2021; St John et al., 2019). A reduced metabolic model of A. niger was constructed by integrating the central metabolism of iJB1325 (Brandl et al., 2018) and beta-alanine pathway for 3-HP production and removing reactions with zero flux at the reference state. The resulting model comprised 172 reactions and 171 metabolites. Proteomics and metabolomics measurements from 51 samples representing 17 distinct strains were used as the observed data. Strain-specific uptake and excretion rates were calculated for measured extracellular metabolites. These rates were further used to calculate internal fluxes with global proteomics as constraints with a proton export objective using E-Flux2 (Kim et al., 2016). At steady state, the metabolic reaction rate ( v ) can be expressed as a function of enzyme concentrations ( e ), internal and external metabolite concentrations ( x and y ) as follows: v = diag v * e e * 1 n + ε x * ⁡ log x x * + ε y * ⁡ log y y * where the asterisks ( * ) denote quantities at a defined reference state, ε x * and ε y * are kinetic parameters called metabolite elasticity matrixes, which respectively quantify the sensitivity of reaction rates to changes in internal and external metabolite concentrations. This lin-log approximation is valid near the reference state and enables efficient parameter estimation (St John et al., 2019).

Posterior distributions of the model parameters were inferred using automatic differentiation variational inference (ADVI), as implemented in the PyMC3 Python library. The model was optimized using the Adagrad optimizer until convergence of the negative evidence lower bound score. Using the resulting parameterized kinetic model, we computed flux control coefficients (FCCs), which quantify the sensitivity of steady-state fluxes to perturbations in enzyme levels at the reference state. FCCs were considered significantly different from zero if their 95% highest posterior density intervals from the posterior distribution do not overlap with zero. Github repositories for the Aspergillus-specific analysis are available at https://github.com/agilebiofoundry/aspergillusq4milestone and for the BMCA at https://github.com/agilebiofoundry/bayesian-metabolic-control-analysis.

Metabolomic and Proteomic analysis was conducted to identify reactions controlling flux in A. niger strains engineered to produce 3-HP. Strains were constructed or selected from a previous studies (Dai et al., 2023) to perturb flux through the beta-alanine 3-HP biosynthetic pathway enabled by overexpression of three heterologous enzymes, TcPAND, BcBAPAT, and EcHPDH (Borodina et al., 2015). Previously identified genetic targets including oahA, uga2, ald6a, ald6b, pyc1, aat1, mct1, and the heterologous beta-alanine pathway genes (Dai et al., 2023; Pomraning et al., 2021) were disrupted and overexpressed to establish a dataset spanning a wide range of 3-HP productivity levels amenable to learn from (Table 1). All strains were cultivated in mRDM for 7 days with extracellular metabolites and biomass collected at days 3, 5, and 7 to provide boundary constraints for flux balance analysis (Figure 1). Cell pellets collected on day 5 were extracted using MPLEx (Nakayasu et al., 2016) for mass-spec based multi-omics assays. Data from 17 strains in triplicate was collected for 55 intracellular metabolites (Kim and Heyman, 2018) and 3,814 proteins using global untargeted methods and 59 proteins selected for quantitatively precise measurement by selected reaction monitoring (SRM) based targeted proteomics (Gao et al., 2020) (Additional file 2).

Intracellular reactions fluxes were predicted using a reduced metabolic model (172 reactions, 171 metabolites) adapted from a genome-scale metabolic model of A. niger (Brandl et al., 2018) that contained only non-zero fluxes at the reference state. The spent media, time, and biomass at collection were used to construct a simple exponential growth model of the organism and estimate strain-specific uptake and excretion rates for key measured extracellular metabolites. These rates were further used to calculate internal fluxes with global proteomics as constraints using E-Flux2 (Kim et al., 2016).

We next employed the BMCA methodology developed to predict how cellular kinetics respond to genetic changes (McNaughton et al., 2021; St John et al., 2019). In BMCA, a low-fidelity kinetic model of microbial metabolism is constructed leveraging linear-logarithmic kinetics. With known kinetic parameters, a kinetic model enables the expected steady-state internal metabolite concentrations and metabolic fluxes to be estimated as a function of enzyme expression and media conditions. With measurements of both the input variables (extracellular metabolite concentrations and enzyme expression) and the output variables (steady-state fluxes and internal metabolite concentrations), posterior distributions in the kinetic parameters that are consistent with the observed data can then be estimated.

Due to the size of the kinetic model considered, posterior distributions in kinetic parameters as a function of the observed data was estimated using automatic differentiation variational inference as implemented in the PyMC3 Python library. The model was optimized until convergence of the evidence lower bound score using the Adagrad optimizer. The posterior predictive distribution (PPD) of the model shows the ability of the model to reproduce the measured experimental data across different strains. The PPD of the fitted model closely reproduces the measured targeted proteomics, metabolite concentrations, and global proteomics-based Eflux2-predicted intracellular fluxes within the unclipped shaded region (Figure 2A). Outside this region, predicted metabolomics, proteins and fluxes were based on clipped measurements, hence the horizontal cluster of sample points.

With a kinetic model and estimated probability distributions in kinetic parameters, we then conducted the Metabolic Control Analysis portion of the BMCA framework. Here, the uncertainty in the estimated kinetic parameters is propagated to the metabolic design strategies suggested by Metabolic Control Analysis. In Figure 2B, we show a subset of the highest posterior density regions of flux control coefficients (FCCs) on 3-HP export calculated from the posterior distribution. A positive FCC indicates that an increase in the corresponding enzyme concentration will increase 3-HP flux, while a negative FCC indicates that a decrease in enzyme concentration will increase 3-HP flux, thus FCCs capture the systems-level regulation of changing enzyme concentration on steady-state metabolic flux.

To evaluate use of the BMCA framework as a predictive modeling tool for synthetic biology applications we selected reactions to modify based on their associated FCC in the training data. We predicted that reactions with negative FCCs will increase flux through the 3-HP pathway when downregulated or deleted and reactions with positive FCCs will increase flux through the 3-HP pathway when enzymes associated with the reaction are over-expressed (Figure 2). The largest positive flux control coefficients include plasmas membrane nitrate transport (r1086, ProteinID 1189116), 5-formyltetrahydrofolate deformylase (r95, ProteinID 1182700), malate dehydrogenase (r44, ProteinID 1144118), pyruvate carboxylase (r19, ProteinID 1031996), and mitochondrial aspartate aminotransferase (r258 m, ProteinID 1184650). The most negative FCC’s include the putative reaction catalyzed by ALD6 (malonate semialdehyde + CoA + NAD(P)⁺ → Acetyl-CoA + CO₂ + NAD(P)H, ProteinID 1182225) (Dai et al., 2023; Pomraning et al., 2021), mitochondrial alcohol dehydrogenase (r118a, ProteinID 1145368) and succinyl-CoA ligase (r38a, ProteinIDs 1145655 and 1141712).

We sought to delete or overexpress each predicted gene target in the highest performing 3-HP strain evaluated in the training set (strain ABF_008348). This strain has a preexisting ald6a disruption, the most strongly suggested deletion target, which we previously defined as a critical modification in 3-HP producing Aspergilli as well as Rhodosporidium toruloides (Dai et al., 2023; Pomraning et al., 2021; Liu et al., 2023). To evaluate alcohol dehydrogenase we selected a high-quality mitochondrial target, adhD (Protein ID 1145368) and assessed both subunits of the succinyl-CoA ligase enzyme. All the single gene overexpression targets were selected for evaluation except plasma membrane transport of nitrate as we had previously established that increasing the concentration of extracellular nitrogen promotes significantly higher yield and titer of 3-HP (Dai et al., 2023). Overexpression of mitochondrial aspartate aminotransferase AAT1 unexpectedly dramatically reduced 3-HP productivity, however, overexpression of a cytosolic variant (AAT2; ProteinID 1176455) improved both titer and yield of 3-HP. In total, six reactions with significant FCCs were evaluated for their impact on flux toward 3-HP. Predicted directional perturbation of five of the six reactions significantly improved titer and rate of 3-HP production and two significantly improved yield (Figure 3). Notably, disruption of the nonintuitive target succinyl-CoA ligase resulted in a 39% increase in titer (14.7–20.4 g/L 3-HP) and 29% increase in yield (0.24–0.31 g 3-HP/g glucose), highlighting the ability to predict meaningful genetic targets in unexpected areas of metabolism by integrating systems-level proteome and metabolome data with metabolic modeling and flux-omics.