Cyanobacteria are the only known prokaryotes capable of oxygenic photosynthesis (Mulkidjanian et al., 2006; Lau et al., 2015). The gram-negative bacteria exhibit a large ecological variety as well as a broad morphological diversity (Bennett and Bogorad, 1973; Rippka et al., 1979; Schirrmeister et al., 2013). Their physiological diversity makes them promising biological chassis for the synthesis of a variety of natural products, including bioactive metabolites like cytotoxins and potential pharmaceutical lead compounds, food supplements, animal feed, pigments, as well as biofuels (Pulz and Gross, 2004; Hays and Ducat, 2015; Jain et al., 2017; Hudson et al., 2021; Barone et al., 2023). Their ability to convert sunlight and atmospheric CO₂ directly into valuable organic compounds could make the chemical and pharmaceutical industry more sustainable and therefore mitigate climate change if high production yields are achieved (Choi et al., 2020; Posten and Schaub, 2009; Oliver and Atsumi, 2015; Oliver et al., 2016).

Metabolic engineering tools for the production of desired compounds are particularly well established for laboratory model strains like Synechocystis sp. PCC 6803 (hereafter Synechocystis) (Knoop and Steuer, 2015; Ramey et al., 2015). The unicellular organism was the first entirely sequenced cyanobacterium (Kaneko et al., 1996) and is one of the best characterized model organisms regarding cyanobacterial biosynthesis (Pils and Schmetterer, 2001; Leplat et al., 2013). Synechocystis is easy and inexpensive to cultivate, and is genetically modifiable with high success rates and predictability (Rippka et al., 1979; Berla et al., 2013).

One promising class of compounds to be produced in cyanobacteria are terpenoids, a large heterogeneous group of naturally occurring organic carbon compounds with over 80,000 known structures (Karunanithi and Zerbe, 2019), having applications in nutrition, medicine and chemistry but also as potential biofuels. Triterpenoids are a group of secondary metabolites, which are composed of 6 isoprene units (Moss et al., 1995), existing in a huge variety of structures with nearly 200 distinct triterpene skeletons, all deriving from the precursor squalene (Xu et al., 2004; Connolly and Hill, 2010). Squalene is typically extracted from shark liver oil, but this method poses serious ecological risks and is not sufficient to sustainably meet increasing demands (Gohil et al., 2019; O’Hagan et al., 2021; Mendes et al., 2022). The diverse applications of squalene include its use as an ingredient in cosmetic products (Gohil et al., 2019), as an antioxidant (Kohno et al., 1995) and an emulsion adjuvant in vaccines (O’Hagan et al., 2021). It has recently been used in several COVID-19 vaccines (Liang et al., 2020; Ho et al., 2021), introducing a surge in demand for this terpenoid. Other reported properties of squalene include tumor-suppressing (Smith et al., 1998; Yang et al., 2014), immunity improving (Ronco and De Stéfani, 2013), cholesterol-lowering (He et al., 2003), as well as antibacterial and antifungal effects (Katabami et al., 2015). Squalene has attracted attention as a feasible source of biofuels (Hellier et al., 2013) as well, if it could be produced sustainably and in large quantities.

In most plants, algae and prokaryotes, terpenoids can be synthesized via the methyl-erythritol-4-phosphate (MEP) pathway, also called the non-mevalonate pathway (Okada and Hase, 2005; Sawai and Saito, 2011). While the mevalonate pathway is present in most eukaryotic cells, the MEP pathway was acquired through endosymbiotic or horizontal gene transfer in plastid-bearing organisms. As the progenitors of plastids, cyanobacteria can serve as both a model organism for chloroplastic terpenoid synthesis through the MEP pathway and are promising production hosts for plant terpenoids (Hemmerlin et al., 2006; Loeschcke et al., 2017). The MEP pathway produces isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP), which are the universal precursors for terpenoid synthesis. In Synechocystis, a single gene, crtE, is responsible for the elongation of terpene precursors towards geranyl pyrophosphate (GPP), farnesyl pyrophosphate (FPP) and geranylgeranyl pyrophosphate (GGPP) in consecutive condensation reactions. The enzyme squalene synthase (Sqs) catalyzes the condensation of two molecules of FPP to presqualene diphosphate (PSPP), which is then converted into squalene via reduction by NADPH. In Synechocystis, squalene is then cyclized by squalene hopene cyclase (Shc) to hopene (Englund et al., 2014). Englund and colleagues (Englund et al., 2014) used a modified strain of Synechocystis to produce squalene by inactivating the shc gene. By preventing the generation of hopene, squalene is accumulated in this mutant. We previously constructed a markerless deletion of shc in Synechocystis to minimize the number of antibiotic resistances carried by each strain and enable further engineering (Dietsch et al., 2021). The availability of inducible promoter systems, such as the rhamnose, anhydrotetracycline or copper inducible promoters in Synechocystis allows fine-tuning of gene expression levels, making improvement of metabolic pathways possible by identifying optimal expression levels for each involved gene (Behle et al., 2020).

A straightforward approach to achieve higher yields of desired products is the overexpression of certain genes to increase the flux towards these metabolites, with many strategies already reported for increasing heterologous terpenoid production (Klaus et al., 2022). Despite some strategies proving successful, the regulation and bottlenecks of the MEP pathway are still not entirely understood (Klaus et al., 2022). Most initial pathway modification approaches result in relatively low product yield, and optimization is often dependent on heuristic techniques (Choi et al., 2017; Steuer et al., 2012). While the identification of genes to be modified is an essential step in metabolic engineering for strain improvement toward the enhanced production of desired bioproducts, it is still difficult to decide which genes to insert or modify, due to the vast number of possibilities of potential targets and the consideration of complex regulation of metabolic networks. In order to rationally identify and overcome bottlenecks for the improvement of strain designs, the use of in silico models has gained increasing significance over the past years (King et al., 2015; Broddrick et al., 2016; Lin et al., 2017; Hendry et al., 2020). Genome-scale metabolic modeling requires only stoichiometric information and no kinetic parameters, which are usually not available even for small reaction networks. In contrast, stoichiometric information is readily and reliably available for a large number of annotated genes (Raman and Chandra, 2009). Additionally, the estimation of stoichiometric yield is computationally not expensive and feasible even for large models involving hundreds of reactions (Knoop and Steuer, 2015). Constraint-based reconstruction and analysis (COBRA) utilizes these models to calculate flux distributions under certain environmental or internal conditions (Schellenberger et al., 2011). COBRA methods are mainly used for the prediction of maximum theoretical yield for native and non-native pathways (Shen and Liao, 2013; King et al., 2015), the effect of gene deletions on biomass or other target compounds (e.g. OptKnock) (Burgard et al., 2003) as well as the identification of bottlenecks and potential targets for up- and downregulation in order to increase product yield (e.g. flux scanning based on enforced objective flux (FSEOF) or minimization of metabolic adjustment (MOMA)) (Segrè et al., 2002; Choi et al., 2010). Previous studies (Choi et al., 2010; Englund et al., 2018; Park et al., 2018), aiming to increase terpenoid production, have shown that constraint-based flux balance analysis (FBA) can be a helpful tool to not only understand and analyze metabolic pathways, but to identify bottlenecks and narrow down the options of potential amplification or knock-down targets for increased product yield, without having a negative influence on the growth rate. Englund et al. successfully employed genome-scale metabolic flux analysis to identify amplification targets that increased isoprene production in Synechocystis (Englund et al., 2018).

In our work, we used an algorithm called FSEOF to screen for potential overexpression targets increasing the terpenoid concentration in Synechocystis (Choi et al., 2010). In silico analysis predicted an increased flux through reactions of the MEP pathway, terpenoid synthesis, the light-dependent reactions of photosynthesis as well as the central carbon metabolism. Additionally, metabolic modeling proposed elevated requirements for the efficient supply, balance and regeneration of cofactors. Twelve of the identified amplification targets were tested in vivo to assess the biological relevance of the model, all of which positively impacted squalene production. The overexpression of sqs led to the strongest increase, with a yield of 13.72 mg l^-1, the highest production titer reported for squalene in Synechocystis to date.

In the pursuit of identifying potential amplification targets, improving cyanobacterial squalene production upon upregulation, we used a genome scale constraint-based metabolic model. Twelve of the identified overexpression targets were experimentally confirmed.

The analysis suggested an up-regulation of fluxes through the MEP pathway and terpenoid synthesis. Since terpenoids are exclusively synthesized via the MEP pathway in Synechocystis, this is very intuitive and shows the robustness of our in silico analysis pipeline. Flux through the light-dependent reactions of photosynthesis as well as the CBB cycle were proposed to increase, to meet the demand of fixed carbon and produce the cofactor NADPH. In fact, the FNR reaction showed an increased flux, while the ATPase reaction is proposed to be downregulated, leading to an overall decreased ATP/NADPH ratio. Previous in silico studies suggested lower ATP/NADPH ratio requirements for many biofuels as well (Erdrich et al., 2014; Shabestary and Hudson, 2016; Englund et al., 2018). Artificially creating imbalances in this ratio could increase product formation. This can be achieved via different approaches, e.g., by blocking cyclic and other alternative electron flows. As a result, ATP and NADPH are forced to be synthesized exclusively via linear electron flow, whose ATP/NADPH ratio is below the ratio required for biomass synthesis. Knocking out ATP producing reactions or introducing ATP futile cycles or other waste reactions could force the organism to use terpenoid production as a sink for excess reduction equivalents (Erdrich et al., 2014). The carboxylation and oxygenation reaction of RuBisCO were coupled in our model, so the flux through the photorespiratory reactions was forced to increase as well. The model suggested using the glycerate photorespiratory bypass. Nonetheless, since it is an undesired reaction to occur, we did not consider it as an amplification target. On the contrary, Zhou et al. could show that impairing photorespiration leads to a redirection of excess energy and isoprene production could be doubled (Zhou et al., 2021). Additionally, cytidine triphosphate (CTP), a cofactor in the MEP pathway, has to be regenerated by multiple phosphorylations of cytidinmonophosphat (CMP) and the model proposed it might become rate limiting upon high squalene flux. Furthermore, the conversion of diphosphate to monophosphate was predicted to be upregulated because of the increased levels of diphosphate released by IspD, CrtE and Sqs. These results are in accordance with previous in silico studies aiming to increase terpenoid production in cyanobacteria, although different models as well as different algorithms were utilized (Lin et al., 2017; Englund et al., 2018). The TCA cycle was indicated to decrease in flux, since it is a sink for carbon, which also was found by previous studies in cyanobacteria (Englund et al., 2018). This finding is in contrast to studies in E. coli (Choi et al., 2010), where the TCA cycle is called to be upregulated. This difference can be explained by the fact that heterotrophs like E. coli have to regenerate their cofactors via the TCA cycle and autotrophs, in contrast, are able to generate cofactors via the light reactions of photosynthesis. The metabolic model suggested a decrease in photosynthetic pigments as well, since increased squalene production diverts carbon flux away from pigments and they compete for the precursor GPP (Lin et al., 2017).

Flux variability analysis found 31 fluxes to be consistent with the results of the FSEOF analysis and increase upon enforced squalene flux. Among the remaining eight fluxes, three showed a pattern where the maximum flux decreased, while the minimum flux increased, and were thus not considered as amplification targets (Choi et al., 2010). The reaction catalyzed by Idi was the only one to display a sign change, where the maximum flux is positive and the minimum flux is negative, implicating the direction of the metabolic flux is unknown and cannot be determined by flux variability analysis (Flowers et al., 2018). Incorporating kinetic parameters as constraints would help to increase the accuracy of the prediction (Moulin et al., 2021). Since in either direction, the flux through Idi positively correlates with squalene production, we chose to verify the result experimentally. Four reactions were completely unbound, thus providing no further information. Previous studies (Choi et al., 2010) showed that the overexpression of unbound targets also had a positive impact on product synthesis.

To confirm the predictions made by the FSEOF analysis, 11 genes were experimentally overexpressed in Synechocystis Δshc. All overexpression strains showed elevated squalene levels, but the degrees to which squalene concentrations were increased varied widely. Overall, growth of the overexpression strains was slightly reduced compared to the control ( Figure S5 (SI)), possibly as a result of the additional antibiotic present in the growth medium. The differences between the overexpression of genes in the MEP pathway are of particular interest, as regulations and feedback mechanisms in the pathway are not entirely known. Since protein tags can affect both activity and stability of enzymes, we chose to use the native protein sequences for overexpression. Since no antibodies are reported for these native enzymes, we were not able to quantify their protein concentrations. However, we confirmed the functionality of our expression system via qRT-PCR for the dxs and sqs overexpression strains which both showed increased transcript levels upon induction with rhamnose ( Figure S4 (SI)). Quantitative differences in RNA abundance between genes are unavoidable due to differing mRNA lengths and stabilities, but it can be safely assumed that induction of expression through addition of rhamnose yields higher gene expression levels than the control in all cases. Since no protein amounts could be quantified in this study, no conclusions can be drawn regarding the metabolic efficiency of each overexpression relative to the expression strength. Instead, changes in the final metabolite concentrations are caused by the combined effect of the strength of the overexpression and the catalytic activity of the protein.

Dxs is reported to be the rate limiting enzyme of the MEP pathway in many previous studies, but did not prove to be a particularly beneficial overexpression target in this study. This may be attributed to the native Sqs activity not providing a strong enough carbon sink downstream of the MEP pathway, leading to the accumulation of intermediates, such as IPP and DMAPP, which are reported to act as inhibitors to the Dxs enzyme (Banerjee et al., 2013; Álvarez-Vasquez et al., 2021; Di et al., 2022). Kudoh et al., 2017 showed that overexpression of dxs led to aggregation of the inactivated protein via allosteric inhibition by IPP and DMAPP, ultimately diminishing the impact of the overexpression on the protein level (Kudoh et al., 2017). Protein inactivation may also explain the small effect of dxs overexpression on metabolites in this study.

Other studies reported MEcPP to accumulate upon overexpression of dxs, so a dual overexpression with ispG may show a more positive effect (Gao et al., 2016; Volke et al., 2019). The overexpression of ispG alone did not lead to increased production however, the dependence of IspG on reduced ferredoxin units may be a limiting factor to the conversion of MEcPP to HMBPP (Wang and Oldfield, 2014). If the concentrations of MEcPP were too low for IspG to be limiting, the binding of ferredoxin to the additional enzyme may be the cause of the reduction in photosynthetic pigments, with the overexpression of ispG leading to the lowest chlorophyll concentrations among the genes of the MEP pathway.

In contrast to most overexpressions, ispE and ispH showed a positive effect on squalene production even without inducer present, but not strongly increasing upon induction. This relationship suggests that these enzymatic steps represent a metabolite bottleneck, which can be relieved through a slight overexpression and then yields diminishing returns upon stronger overexpression.

The overexpression of idi led to an increase in squalene while decreasing pigments, which is in accordance with the predictions made by the constraint-based model. IspH favors production of IPP from HMBPP over DMAPP, while Idi favors the isomerization of IPP towards DMAPP, leading to an IPP : DMAPP ratio of 3:1 in vivo, the optimal ratio for production of GGPP (Chaves et al., 2016; Volke et al., 2019). The overexpression of idi likely shifted the IPP : DMAPP ratio in favor of DMAPP.

Overexpressing sqs led to the strongest increase in squalene and reduction in carotenoids. Sqs competes with the native GGPP synthase crtE for the intermediate FPP, with the overexpression leading to a shift in favor of squalene, away from GGPP, the precursor for carotenoids and phytol. Since there was only little change in growth over three days ( Figure S5 (SI)), the reduced pigmentation did not seem to have a significant negative effect on the cell, but shifting the balance further towards squalene might lead to decreased photosynthetic performance. Overexpression of ispH and ispE on the other hand increased squalene concentrations as well as the total amount of measured terpenoids, defined as the sum of squalene, carotenoid and the phytol chain of chlorophyll by up to 18%. Combining the overexpressions of sqs with ispH and ispE, may be a promising strategy moving forward, as increased total flux through the MEP pathway can both compensate for the reduced pigmentation caused by sqs overexpression and increase squalene titers.

In conclusion, FSEOF allowed us to choose biologically relevant amplification targets computationally, all of which had a positive effect on squalene synthesis upon experimental validation. Considering squalene is synthesized via the linear MEP pathway, most of the identified targets are rather intuitive. Since we were able to confirm all selected targets, we suggest a validation of the non-intuitive targets outside the MEP pathway, such as Ppa and FNR for further studies. These two targets have previously been predicted and tested in Synechocystis (Englund et al., 2018). Our findings propose that constraint-based metabolic models could aid in the selection of targets improving the production of desired metabolites. This could be of particular interest for the prediction of combinatorial interventions. However, classic FBA does not account for regulatory mechanisms like feedback inhibition or potential metabolite toxicity (Knoop and Steuer, 2015). It could be helpful for future studies to test FSEOF in combination with dynamic extensions of FBA (Mahadevan et al., 2002) and hybrid kinetic and constraint-based models (Shameer et al., 2022). The incorporation of metabolite concentrations, flux rates or kinetic parameters could drastically improve the precision and reliability of the results (Mahadevan et al., 2002; Moulin et al., 2021; Shameer et al., 2022). For future studies we suggest the combinatorial expression of the identified amplification targets, especially with the native sqs, since it seems to be the rate limiting enzyme in the present study, as well as the combination of amplification and knock-down targets. Our experiments are in accordance with the central idea of synthetic biology, where experimental designs are determined by metabolic modeling and experimental results can feed back data into models to increase their accuracy, leading to deterministic, steady improvement. The overexpression of the native sqs gene of Synechocystis proved to be the most successful strategy for squalene production to date, with the strain reaching a higher production titer than heterologous sqs expression in Synechocystis (Pattanaik et al., 2020).

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

AG: Conceptualization, Investigation, Methodology, Computational work, Writing – Original draft, review & editing, Data Visualization. AN: Conceptualization, Investigation, Experimental work, Writing – Original draft, review & editing, Data Visualization. MD: Conceptualization, Methodology, Experimental work, Writing – review & editing. TP: Computational work, Writing -Data visualization, review. DM: Experimental work, Writing – Review & Editing. ST: Experimental work, Writing – Review & Editing. PW: Methodology, Data Curation. IA: Supervision, Writing – review & editing. All authors contributed to the article and approved the submitted version.