Crassulacean acid metabolism is a photosynthetic adaptation syndrome of carbon concentration that has evolved multiple times in plants (Crayn et al., 2004; Winter and Smith, 2012). It consists in the nocturnal fixation of CO₂ into malate, which is stored in the vacuole and decarboxylated the next day providing CO₂ for Rubisco (Dodd et al., 2002). It is commonly associated with conditions of water scarcity in which plants close their stomata during the day to avoid water loss. Keeping stomata closed during the day while concentrating CO₂ at night permits a high water use efficiency, allowing plants to remain productive even under dry conditions. A great interest in CAM crops has recently developed, given the potential of extending their use in semi-arid and degraded lands (Liu et al., 2018). In natural ecosystems, CAM plants are less affected by increased nighttime temperatures and rainfall variability than C4 grasses (Huang et al., 2020) and are increasing their cover in many arid regions (Yu and D’Odorico, 2015). Global climate models predict drier regimes and decrease in soil moisture in large areas of the world (Dai, 2013). In many of today’s agricultural lands, CAM crops could be the only crops able to grow in a fairly close future. Understanding CAM evolution is especially relevant for CAM engineering into C3 crops. The multiple occurrences of CAM in numerous lineages across the plant kingdom and different levels of CAM expression suggest that the evolutionary path from C3 to CAM is relatively straightforward (Schiller and Bräutigam, 2021). Finding the key events that led to CAM emergence can potentially pinpoint the key transformations necessary for engineering CAM into C3 plants.

Mathematical models on CAM metabolism provide valuable information. Two types of models have been used to study CAM: ordinary differential equations (ODE) and flux balance analysis (FBA) models. ODE-based models can accurately describe metabolite concentrations and their change over time at the cost of requiring quantitative knowledge on the parameters and enzyme activities of the pathways being modeled. On the other hand, FBA-based models require no quantitative a priori knowledge. They can be used to study steady state flux distributions but provide no information on concentrations or transient dynamics.

ODE-based models have determined the components and regulatory interactions sufficient to bring about CAM behavior, while FBA models have shown that many CAM defining features —such as nocturnal malate accumulation—are the result of flux optimization upon limited CO₂. Here, we review the findings of both groups of models. Although most of the works were not conceived with evolution in mind, several of their findings have direct evolutionary implications, in particular for the ongoing discussion on the C3-to-CAM evolutionary continuum (Cushman and Bohnert, 1999; Dodd et al., 2002; Silvera et al., 2010; Bräutigam et al., 2017).

Crassulacean acid metabolism is divided into four archetypical phases. Phase I occurs at night when the stomata are open and starch reserves are used to synthesize phosphoenolpyruvate (PEP). High amounts of PEP are needed to fix CO₂ overnight, a reaction catalyzed by PEP carboxylase (PEPC). PEP carboxylase is reversibly phosphorylated during the night to reduce malate feedback inhibition (Carter et al., 1991). Oxaloacetate (OAA) is the product of PEPC and is quickly converted by malate dehydrogenase into malate (Cuevas and Podestá, 2000) which is transported into the vacuole. A vacuolar H⁺-ATPase moves large quantities of H⁺ into the vacuole generating a proton gradient that energizes malate transport (Cheffings et al., 1997). Phase II occurs during the first hours of dawn. In this phase, the transition from PEPC-mediated to Rubisco-mediated CO₂ fixation occurs, and for a period of time CO₂ is fixed by both carboxylases (Lüttge, 2015). Phase III takes place during the day while stomata are closed. During this phase, CO₂ is freed from malate and re-fixed by Rubisco. For this to happen, malic acid is re-mobilized by passive efflux and decarboxylated. There are two alternative pathways for malate decarboxylation, depending on the species. One starts with the release of CO₂ by the conversion of malate to pyruvate by NAD(P)-ME, followed by the regeneration of PEP by PPDK (Kondo et al., 2000). In a second pathway, NAD(P)-MDH converts malate to OAA, followed by PEPCK, regenerating PEP, and releasing CO₂ (Cushman and Bohnert, 1999). Finally, phase IV takes place late in the afternoon once malate reserves have been depleted. At this time, stomata open again and allow for the uptake of CO₂, which is processed directly by Rubisco (Dodd et al., 2002).

The occurrence of CAM four phases varies from species to species and depends on the developmental stage or the severity of environmental stress. Young tissues of all CAM plants perform C3 photosynthesis (Winter, 2019). Upon drought stress, phases II and IV can be lost or reduced in constitutive CAM plants (Dodd et al., 2002). The proportion of PEPC-fixed CO₂ also increases with drought in species with a C3-CAM intermediate metabolism, while facultative CAM species can revert completely to C3 after rewatering (Winter and Holtum, 2014). In addition, several CAM modes are possible; including CAM cycling, consisting in respiratory CO₂ re-fixation, but with stomatal opening during the day; archetypical CAM; and CAM idling, which happens upon severe drought, in which only respiratory CO₂ is re-fixed and stomata remain closed day and night (Borland et al., 2011).

The functioning of CAM metabolism results in the diurnal oscillation of malate content and CO₂ uptake. As the oscillations continue despite leaving the plants in continuous light, or darkness and CO₂-free air (Wilkins, 1959, 1960; Nuernbergk, 1960), it became clear that understanding these dynamics would require a mathematical framework. Nungesser et al. (1984) translated a basic scheme of reactions known to be relevant for CAM metabolism into a set of ordinary differential equations (ODE). This so-called skeleton model would become the basis for a series of models that would continue for two more decades, dealing with the perturbation of the rhythmicity as well as with finding the dynamically essential components (Table 1).

Nungesser et al. (1984) skeleton model included six regulatory steps: inhibition of PEPCase by malate, activation of this last reaction by Glc6P, inhibition of starch consumption by PEP, activation of photosynthesis by light intensity, inhibition of CO₂ uptake by increased internal concentration of CO₂ (C_i), and activation of malate transport from the vacuole to the cytosol by light (Figure 1). The predicted relative content of Glc6P, PEP, malic acid, starch, C_i, and net CO₂ exchange along the 24 h cycle were remarkably in agreement with experimental data. This shows that the enzymes and metabolite pools considered in the model are sufficient to explain the basic features of CAM.

Nungesser et al. (1984) assumed that the transport of malate between the vacuole and the cytoplasm was stimulated by light. However, oscillations in malate content under continuous light challenged this idea. To solve this issue, they added an oscillating element—a hysteresis switch—controlling malate transport through the tonoplast that would be used in the future by other authors. Lüttge and Beck (1992) refined the model by adjusting several parameters to stabilize oscillations under prolonged runs, kept the hysteresis switch, and eliminated the direct influence of light on malate transport. The changes allowed the simulations to reproduce the endogenous malate oscillations under continuous light.

Three subsequent refinements improved model fitting to the observed effect of irradiance and temperature on oscillations (Blasius et al., 1997, 1998, 1999). The original six pools of metabolites were reduced to only three: cytoplasmic and vacuolar malate pools, and C_i. This model assumes an active influx filling the vacuole with malate during the night. At a critical concentration, the hysteresis switch kicks in and vacuolar malate is rapidly expelled. C_i is the carbon source for malate production, while malate consumption and CO₂ assimilation connect the oscillatory behavior with the rest of metabolism. Blasius et al. (1999) improved the mechanistic details that generated the hysteresis switch by considering membrane properties and simulating tonoplast permeability as a function of membrane order.

The precise malate transporter responsible for the hysteresis switch remains to be identified. Meyer et al. (2011) found that malate channel ALMT6 shifts from inward to outward rectification in Arabidopsis guard cells at a low vacuolar pH. It was recently suggested that a similar mechanism could operate in CAM plants (Ceusters et al., 2021). Thus, it might be a drop in pH and not a critical vacuolar turgor that triggers the influx-to-efflux shift to begin malate decarboxylation. Transcriptomics data from the CAM species Agave americana suggest the malate/citrate antiporter tDT might also play a role in overnight malate accumulation (Yin et al., 2018).

The fact that skeleton models are sufficient to reproduce CAM dynamics suggests that CAM can originate from the emergence of such metabolic circuits during evolution. After a quick literature survey based on the model by Nungesser et al. (1984), we found evidence that not only enzymes but also regulatory interactions are already in place in C3 plants (Figure 1). These regulatory interactions play roles in different primary metabolism pathways and date back to different stages in the evolution of life. PEP inhibition of glycolytic enzymes is a feedback regulation mechanism that is shared with heterotrophic bacteria and animals (Plaxton, 1996; Ogawa et al., 2007; Grüning et al., 2014). Malate inhibition of PEPC has been documented in unicellular green algae and is crucial for coordinating carbon metabolism with nitrogen assimilation (Peak and Peak, 1981; Schuller et al., 1990). The activation of photosynthesis by light and the subsequent carbon fixation by Rubisco is common for all photosynthetic organisms (McKay et al., 1991; Portis et al., 2008; Tsai et al., 2015). Stomatal closure induction by high C_i is known to occur in early vascular plants and is thought to have been essential for the colonization of land by plants (Chater et al., 2011). Glc6P activation of PEPC is found in C3, C4, and CAM plants (Chollet et al., 1996; Blonde and Plaxton, 2003). In C3 plants, this regulation plays a role in gluconeogenesis by regulating the rate of starch mobilization with the production of dicarboxylic acids (Law and Plaxton, 1995). Finally, there is evidence that the autonomous vacuolar oscillator is also present in C3 plants. As malate can be toxic for the cell in large quantities and is essential for a number of biochemical pathways, its amount in the cytosol is tightly regulated by storing excess malate in the vacuole (Fernie and Martinoia, 2009). Vacuoles of C3 plants switch from malate influx to efflux depending on a critical concentration of cytosolic malate (Gout et al., 1993). As in CAM plants, malate in the vacuole of C3 plants fluctuates diurnally (Gerhardt et al., 1987).

The accumulated information on biochemical reactions together with the availability of genomic data made the first FBA genome-wide metabolic models possible (Schilling and Palsson, 1998; Schilling et al., 1999; Edwards and Palsson, 2000). In plants, the first FBA whole-genome (FBA-WG) metabolic reconstruction was done for Arabidopsis (Poolman et al., 2009), followed by maize (Saha et al., 2011). Later, Cheung et al. (2013) incorporated transport and maintenance costs into the Arabidopsis whole-genome FBA model, allowing for a more realistic prediction of flux distribution in central metabolism. Building on this improvement, Cheung et al. (2014) incorporated the interaction between diurnal and nocturnal metabolism (Table 2). Essentially, this was done by creating two alternative models corresponding to each phase. The models were coupled by adding “transport” reactions between them that would pass on accumulated metabolites. A specified amount of sucrose and amino acids was set to be exported to the phloem, thus requiring the cell to be productive. The flux distribution of both models was simultaneously optimized to decrease the overall amount of material flowing through metabolism. Since the study addressed diel metabolic changes, it was natural to use the model to compare C3 and CAM scenarios. This was done by setting CO₂ exchange to zero during the night for CAM. Remarkably, the simulations are in agreement with knowledge on plant metabolism. The C3 simulation predicted citrate to accumulate during the night and to be used during the day for amino acid biosynthesis. For CAM, the nocturnal accumulation of malate was predicted with PEPC as the carboxylating enzyme.

Shameer et al. (2018) simplified Cheung’s FBA-WG model (Cheung et al., 2014) to a core stoichiometric model of central plant metabolism and considered explicitly organellar pH and the charged state of metabolites. The modifications allowed the improved model to match experimental measurements of starch and malate and revealed that CAM metabolism is more productive—higher phloem export rates—than C3, compensating the higher energetic costs of CAM metabolism. These results are highly valuable as a long-sought goal of biotechnologists has been the implementation of CAM into C3 crops by means of genetic engineering.

Shameer et al. (2018) model was the basis for two further works (Töpfer et al., 2020; Tay et al., 2021). Töpfer et al. (2020) extended the model by allowing for differences in vacuolar storage capacity, making gas exchange dependent on temperature and relative humidity, and considering the hourly progression of metabolism during a day. These last two additions allowed for CO₂ uptake variations that resulted from a compromise between optimizing CO₂, metabolic demand, and minimizing water loss. Interestingly, under a C3-like vacuole size restriction, flux optimization leads to reduced CO₂ uptake during the hottest hours of the day—equivalent to midday depression of photosynthesis. A similar simulation under a CAM-like vacuole size restriction increased CO₂ uptake at the end of the night together with carboxylic acid accumulation.

Tay et al. (2021) used the Shameer et al. (2018) model to study the C3-to-CAM continuum. They simulated CAM idling and CAM cycling—both considered as weak CAM—by impeding CO₂ uptake during the night and during the whole cycle, respectively. The modified models predicted flux distribution intermediate between C3 and CAM metabolism in terms of malate and starch accumulation. In a second analysis, the authors performed a series of simulations decreasing CO₂ uptake gradually. They observed a steady shift from C3 to CAM photosynthesis with a respective increase in flux through the starch/malate cycle.

Adaptive evolution reshapes metabolism to maximize fitness (Long and Antoniewicz, 2018). If we accept this, given adequately identified optimality criteria, mathematical models can be useful to understand and predict the course of evolution (Bordbar et al., 2014; de Visser and Krug, 2014). FBA models maximize fitness theoretically by translating a selective pressure into an objective function and have provided remarkably accurate predictions (Varma and Palsson, 1994; Harcombe et al., 2013). The set of FBA models studying the metabolic interactions between day and night suggests that simulating longer periods of stomatal closure are enough to trigger a flux distribution typical for CAM metabolism. Remarkably, CAM-like fluxes emerge upon constraining CO₂ uptake on models derived from a C3 whole-genome metabolic reconstruction. Regardless of whether stomatal closure happens during the day or the night and regardless if the stomatal closure is full or partial, flux distributions resembling CAM metabolism are evident. Therefore, assuming that evolution selects for efficient responses to selective pressures (Harcombe et al., 2013), the models provide formal evidence for limited CO₂ as the main driver of CAM evolution, a hypothesis long proposed and widely accepted (Lüttge, 2002).

Reduced CO₂ uptake can result from limited water availability, as CAM-like fluxes appear when CO₂ uptake is compromised by maximizing water saving (Töpfer et al., 2020). At low air humidity or when roots experience water deficit, different signals such as jasmonates and ABA induce stomatal closure (Bauer et al., 2013; De-Ollas et al., 2018). At hot dry hours, C3 and C4 plants close their stomata and undergo what is called midday depression of photosynthesis (Correia et al., 1990; Hirasawa and Hsiao, 1999). Töpfer et al. (2020) model was able to recapitulate this phenomenon. When water saving was maximized, a prolonged stomatal closure during the day was predicted. It is possible that as plants colonized dry areas or whole regions of the world became arid, the midday depression of photosynthesis widened, narrowing the time period at which CO₂ exchange could occur. Extended hours behind closed stomata could have exerted a significant selective pressure to develop CAM so that plants could remain productive even under limited CO₂ exchange.

A reduction in atmospheric CO₂ has also been proposed to trigger CAM evolution, as CAM radiation events coincide with declining CO₂ in the geological records (Keeley and Rundel, 2003; Silvera et al., 2009). Tay et al. (2021) tested this hypothesis by gradually decreasing the simulated concentration of CO₂ and observed CAM flux distributions as well. A sustained condition of limited CO₂ could have selected for more efficient metabolic fluxes over time that resulted in what we know as CAM. Tay et al. (2021) model provided several predictions for the metabolic transitions during CAM evolution worth testing experimentally by further works.

Although conceived with a different purpose, ODE models addressing CAM make also a significant contribution to the understanding of CAM evolution. Our quick literature survey found that not only the enzymes, but the regulatory interactions in CAM skeleton models are already in place in C3 plants. These regulatory interactions date back to different times in evolution and belong to different pathways in primary metabolism. Borrowing existing reactions is commonplace in the evolution of new pathways since it is faster and energetically less expensive than generating de novo reactions. The Krebs cycle is believed to have arisen from primitive pathways for amino acid synthesis (Meléndez-Hevia et al., 1996). Interestingly, CAM has also been proposed to be derived from an amino acid synthesis pathway (Bräutigam et al., 2017). In C3 plants, the accumulation of organic acids during the night fuels amino acid synthesis the following day (Gauthier et al., 2010). In this pathway, PEPC carboxylase is active during the night and catalyzes the synthesis of PEP resulting from carbohydrate degradation (Nimmo, 2003). In agreement with this, interactions (1), (2), and (3) of CAM skeleton model (Figure 1) can coordinate carbon flux toward amino acid synthesis in C3 plants: Glc6P is a feedforward activator of PEPC, while malate and PEP are feedback inhibitors of PEPC and glycolysis, respectively. Excluding vacuolar storage, this pathway is also found in unicellular green algae, in which PEPC activity is regulated to coordinate carbon metabolism with nitrogen assimilation (Peak and Peak, 1981).

The crucial step that diverts the pathway from amino acid synthesis to a new type of photosynthetic metabolism is malate decarboxylation. Malate decarboxylation connects the pathway regulated by interactions (1), (2), and (3) with carbon re-fixation by Rubisco during the day (4) and the accumulation of CO₂ that induces stomatal closure (5; Figure 1). Either NAD(P)-ME or PEPCK are responsible for the release of CO₂, depending on the pathway. Neither enzyme was originally part of photosynthetic metabolism (Sanwal and Smando, 1969; Boles et al., 1998; Drincovich et al., 2001), but we can presume that mutations that led to the induction of NAD(P)-ME and PEPCK activity during the day could have been evolutionarily advantageous by alleviating CO₂ starvation during extended periods of stomatal closure. Both enzymes have been recruited multiple times during evolution when CO₂ is limiting (Reiskind and Bowes, 1991; Magnin et al., 1997; Casati et al., 2000; Hibberd and Quick, 2002), besides their well-known role in C4 and CAM metabolism (Bovdilova et al., 2019; Lim et al., 2019). For the CAM subtypes using NAD(P)-ME, NAD(P)-ME induction could have been a direct consequence of an initial increase in PEPC activity. In heterotrophic tissues relying primarily on oxidation, PEPC and NAD(P)-ME work in concert to provide reducing power in the form of NAD(P)H (O'Leary et al., 2011). Transgenic C3 plants overexpressing PEPC exhibit an important increase in NAD(P)-ME expression and NAD(P)-ME enzymatic activity, probably triggered by an excess of cytosolic malic acid (Häusler et al., 2001). It has been previously proposed that during CAM evolution NAD(P)-ME and PEPC activities became coordinated to avoid futile CO₂ cycling between both enzymes (Dodd et al., 2002). If high cytosolic malic acid induces NAD(P)-ME, the storage of malate during the night and the switch from malate influx to efflux in the vacuole during the day could theoretically sort out the issue of the overlap of both enzymatic activities, although this possibility has not been tested yet.

In order for CAM skeleton models to run, enough amounts of malate need to accumulate in the vacuole so that decarboxylation builds up a high C_i that induces stomatal closure, and thus, the reverse stomatal patterns appear. Vacuolar size is critical for CAM function (Griffiths et al., 2008). Considering this parameter explicitly, Töpfer et al. (2020) showed that the typical malate storage capacity of C3 plants is not sufficient to allow the emergence of CAM flux distributions. Thus, changes that lead to an increased malate storage are necessary in order to develop a proper CAM pathway, these could include H⁺-ATPase (Ratajczak et al., 1994), malate channels (Pantoja and Smith, 2002), and malate/citrate antiporters (Frei et al., 2018).

So far, the proposed ODE models consider only the control of C_i on stomatal aperture. Although there is substantial evidence for this regulation (Cockburn et al., 1979; Males and Griffiths, 2017), we know now that it is not the only one. In order to undergo the transition from diurnal to nocturnal stomatal aperture, C3 plants need to rewire the circadian control on stomata (Hassidim et al., 2017). Comparative genomic studies have shown that genes involved in stomatal regulation have shifted their peak expression to the evening together with clock or clock-entraining genes (Yang et al., 2017; Yin et al., 2018). Thus, it is also possible that the inverse pattern of stomatal aperture was achieved before attaining a high C_i during the night.

Another theoretical possibility is that the circadian control on stomatal aperture was lost during CAM evolution and stomatal aperture relied mostly on the autonomous vacuolar oscillator (6). This way autonomous oscillations of malate content could drive stomatal closure during the day on their own. More research is needed to understand how autonomous oscillations in malate content and the circadian clock are integrated to control CAM rhythmicity. Currently, it is known that both, impeding the circadian input on CAM and silencing CAM, alter oscillations in CO₂ uptake and the circadian clock itself (Boxall et al., 2017, 2020). Further ODE models integrating the clock with CAM metabolism could not only shed light on CO₂ uptake regulation, but on the overall interplay between carbon metabolism and the clock.

An ongoing discussion is the C3-to-CAM continuum itself (Edwards, 2019; Yang et al., 2019; Schiller and Bräutigam, 2021; Winter and Smith, 2022). On one hand, the remarkable phenotypic plasticity observed in CAM plants generates a wide gradient in δ^13C—an indicator of the amount of carbon fixed by CAM—in several plant lineages (Silvera et al., 2010). On the other hand, there is a clear bimodal distribution showing there are only a few intermediate phenotypes in comparison with C3 and strong CAM phenotypes, which might indicate mechanistic constraints for intermediate metabolism (Borland et al., 2011). Yang et al. (2019) argue that the continuum hypothesis underestimates the number of changes needed to achieve CAM, since available genomics data indicate at least 60 genes are involved in CAM evolution. Winter and Smith (2022) highlight that the change to nocturnal malate accumulation would involve a qualitative leap in CAM evolution since for many C3 species malate accumulates during the day.

The discussion truly appears as another instance of the old saltationism versus gradualism debate (Gould and Eldredge, 1977). However, it has opened the conversation on the pace of CAM evolution and the conditions that would facilitate either a smooth or sharp transition. Edwards (2019) analyzes the subject in terms of evolutionary accessibility and proposes there is smooth transition between C3 and weak CAM phenotypes, but for strong CAM to occur, the development of extreme succulence is necessary. Future modeling efforts on CAM metabolism could provide more evidence to discuss the subject. It would be highly informative to combine genomic information, evolutionary modeling and genome-wide models to investigate possible points of evolutionary contingency. In addition, mathematically controlled comparisons can be used to understand the evolutionary limits of the various metabolic types (Alves and Savageau, 2000; Schwacke and Voit, 2004). Another area of opportunity for CAM models is to aid in the rational implementation of CAM metabolism in C3 crops. Recently, a set of 13 CAM enzymes were overexpressed independently in Arabidopsis (Lim et al., 2019). Remarkably, the lines expressing the carboxylation module resulted in increased stomatal conductance and increased malate accumulation, and most lines from both carboxylation and decarboxylation modules had increased biomass. Feeding FBA models with experimental measurements of these plants can indicate how fluxes may be optimized after the overexpression and provide clues for novel bioengineering strategies.

Computational analysis and modeling of CAM metabolism have not only contributed to the understanding of how CAM basic components integrate for the emergence of an autonomous rhythm, and of how fluxes through central metabolism are reshaped, but may also contribute significantly to identify evolutionary selection processes. The convergence of different FBA models to CAM flux distributions provides formal mathematical evidence for low C_i as the main driver for the evolution of CAM. We believe that the use of these models together with evolutionary models and genomic data will shed more light on CAM evolution and assist rational engineering of CAM into C3 crops.

AB conceived the idea. AB, EM, EV, IM-C, and RA wrote the manuscript. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.