Divergent Transcriptional Architectures Beyond Core CAM Genes in Facultative and Constitutive CAM Species in Tillandsia L
Clara Groot‐Crego, Sarah Saadain, Marylaure De La Harpe, Jaqueline Hess, Michael H. J. Barfuss, Walter Till, Gert Bachmann, Wolfram Weckwerth, Christian Lexer, Ovidiu Paun
TL;DR
This study explores how different species use similar core genes but distinct supporting genes to evolve water-efficient photosynthesis called CAM.
Contribution
The research reveals that CAM evolution involves conserved core genes but divergent transcriptional architectures in closely related species.
Findings
Core CAM enzymes show conserved expression patterns across species.
Facultative CAM species uniquely upregulate PPC2 instead of the canonical PPC1 under water stress.
Transcriptional shifts in stomatal movement and starch metabolism show minimal overlap between species.
Abstract
Crassulacean acid metabolism (CAM) is a water‐efficient photosynthetic strategy involving a coordinated suite of complex traits including metabolic, anatomical and regulatory aspects that shift across the diel cycle. While CAM has evolved repeatedly in land plants, the evolutionary routes enabling this convergence remain elusive. Whereas the same core CAM (de)carboxylation genes are consistently involved, a key question is whether distinct CAM phenotypes also depend on a shared set of auxiliary genes, reflecting a quantitative continuum of expression, or whether they can instead emerge through divergent or redundant peripheral solutions. The bromeliad subgenus Tillandsia, with diverse photosynthetic strategies, offers an ideal system to explore this question. Using physiological and transcriptomic analyses of well‐watered and water‐limited accessions of two closely related species, we…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7
FIGURE 8
FIGURE 9- —Austrian Science Fund10.13039/501100002428
- —Agence Nationale de la Recherche10.13039/501100001665
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsPhotosynthetic Processes and Mechanisms · Plant and Biological Electrophysiology Studies · Light effects on plants
Introduction
1
The repeated evolution of complex traits has been a long‐standing matter of study for evolutionary biologists. It remains unclear how traits which need coordinated changes across multiple organismal layers can evolve frequently and repeatedly across the tree of life. One explanation is that similar phenotypes can arise through different sets of mutations, also referred to as ‘genetic routes’, assuming complex traits are polygenic and the underlying genes show substantial redundancy (Barghi et al. 2020; Goldstein and Holsinger 1992; Láruson et al. 2020; Nowak et al. 1997). Indeed, convergent phenotypes lacking a shared genetic basis have been documented across a broad spectrum of life including beach mice (Steiner et al. 2009), conifers (Yeaman et al. 2016), laboratory‐reared Drosophila (Barghi et al. 2019), cacao plants (Hämälä et al. 2020) and Heliosperma ecotypes (Szukala et al. 2023).
Crassulacean acid metabolism (CAM) represents a prime example of a collection of coordinated, complex yet highly evolvable traits, having evolved independently at least 66 times across 38 plant families (Gilman et al. 2023). CAM plants reconfigure photosynthetic metabolism across the diel cycle by assimilating CO_2_ at night rather than during the day, as C_3_ plants do. Nocturnally fixed CO_2_ is then temporarily converted to malate by phosphoenolpyruvate carboxylase (PEPC) and malate dehydrogenase (MDH). Malate is stored in vacuoles and later decarboxylated during the day to supply CO_2_ for the Calvin–Benson–Bassham cycle. This enables daytime stomatal closure, increasing Rubisco efficiency by raising intracellular CO_2_ levels (Osmond 1978) and drastically reducing water loss through evapotranspiration (Borland et al. 2014). These adaptations are especially advantageous under high temperatures or drought stress, which exacerbate Rubisco oxygenation (Heyduk et al. 2021).
CAM integrates several interdependent traits, including temporally coordinated transcriptional regulation to set the alternate biochemical pathway in motion, and cellular processes like stomatal movement, sugar and malate transport and starch metabolism (Borland et al. 2016; Schiller and Bräutigam 2021). Additionally, certain anatomical features often appear associated with CAM, leading to the hypothesis that they may be a prerequisite for a plant to fully express CAM (Nelson et al. 2005). Examples of such anatomical features are enlarged vacuoles (Heyduk et al. 2016; Males 2018), tighter packed mesophyll cells (Nelson and Sage 2008) and succulence (Gibson 1982; Griffiths et al. 2008), though for the latter, intraspecific variation does not appear correlated with CAM (Martin et al. 2009).
CAM was initially described as having a discrete, bi‐modal distribution of phenotypes, yet detailed physiological studies of C_3_ species with CAM relatives have revealed weak and facultative forms of CAM that occur without major anatomical changes (Heyduk et al. 2021; Neales 1975; Pierce et al. 2002; Silvera et al. 2005; Winter and Holtum 2002). These findings prompted a shift towards conceptualising CAM as a continuum, though debate persists (Bräutigam et al. 2017; Schiller and Bräutigam 2021; Winter and Smith 2022). The CAM continuum concept remains difficult to reconcile with the many discrete CAM categories in the literature (e.g., facultative, weak, strong, constitutive CAM, CAM‐cycling; Winter 2019), which do not clearly align with a continuous phenotype spectrum (Edwards 2023). Additionally, weak/facultative CAM plants can be difficult to detect, as carbon isotope ratios, the most cost‐effective method to measure CAM activity on a large taxonomic scale, lack sensitivity for these phenotypes (Pierce et al. 2002; Winter and Holtum 2002). More detailed measurements, such as titratable acidity or gas exchange under standard and drought conditions, are needed to distinguish weak/facultative CAM from C_3_ species.
The evolutionary trajectory of CAM and the molecular similarity across CAM phenotypes remain unclear. While diel expression shifts of core CAM enzymes such as PEPC appear broadly conserved among CAM plant taxa, the realised CAM phenotype likely involves numerous small‐effect genes stemming from broader metabolic, regulatory and anatomical networks, potentially with a high redundancy. Comparative analyses between different CAM phenotypes can reveal whether CAM evolution follows a constrained, shared evolutionary path or instead involves large genetic redundancy. For example, Heyduk et al. (2018) found both shared and variable expression patterns in CAM‐related genes among weak and strong CAM phenotypes in Agavoideae, but more studies comparing the molecular basis of CAM between closely related taxa are needed to obtain a comprehensive view of genetic redundancy and constraints in CAM evolution.
The subgenus Tillandsia (genus Tillandsia L., Bromeliaceae) is an emerging model for studying repeated and rapid evolution of key innovation traits. As part of one of the fastest diversifying plant lineages (Givnish et al. 2014), the subgenus comprises an adaptive radiation of over 250 species formed within the last 7 million years (Barfuss et al. 2016; Mendoza et al. 2017; Yardeni et al. 2025). Numerous key innovations characterise this radiation, including CAM, which spans a broad phenotypic range from C_3_‐like to strong constitutive CAM (Crayn et al. 2004; De La Harpe et al. 2020).
In this study, we compare a pair of closely related species exhibiting facultative and constitutive CAM to assess the degree of genetic redundancy underlying this suite of interdependent complex traits. CAM is a complex phenotype influenced by multiple key genes and protein interactions, cellular processes and metabolic pathways. If CAM evolution is tightly constrained, both phenotypes should represent positions along a shared evolutionary trajectory, mediated by the same set of genes but expressed to different degrees. Conversely, limited overlap in gene expression, at least in some parts of the network, would suggest a (partly) redundant genetic architecture, consistent with greater evolvability and multiple independent evolutionary routes to CAM evolution.
The species T. vanhyningii Harms (formerly T. ionantha var. vanhyningii L.B.Sm., elevated to species level by Beutelspacher and García‐Martínez (2021)), and T. leiboldiana Schltdl. are closely related members of the subgenus Tillandsia, with a maximum divergence time of 3.2 Mya (Barfuss et al. 2016; Mendoza et al. 2017; Yardeni et al. 2025). Tillandsia vanhyningii is a succulent ‘grey’ species with dense trichome cover, whereas the ‘green’ T. leiboldiana typically forms an impounding tank and is comparatively less succulent. These species display some of the most disparate carbon isotope values in the subgenus (De La Harpe et al. 2020), with T. vanhyningii having δ^13^C values typical of strong, constitutive CAM (−13.9), while T. leiboldiana shows values typical of C_3_ photosynthesis (−31.3). However, recent metabolomic, genomic and transcriptomic analyses suggest that T. leiboldiana may be a facultative CAM species (Groot Crego et al. 2024). Despite minimal night‐time malate accumulation, T. leiboldiana shows CAM‐like expression in certain CAM‐related genes, including PEPC kinase, under well‐watered conditions, indicating a latent CAM capacity potentially activated under drought‐induced stress.
We subjected accessions of both species to standard and water‐limited conditions and measured titratable acidity, gas exchange and transcriptomic profiles across the diel cycle to compare their photosynthetic phenotypes under both watering regimes (Figure 1). Our study characterises their photosynthetic phenotypes and elucidates how water‐withholding modulates CAM expression. We show that, although both species perform CAM under water‐limited conditions, they exhibit substantial divergence in transcriptional modulation of CAM‐related peripheral genes beyond the core CAM genes. These findings illuminate the degree of genetic redundancy underlying CAM and provide insight into the multiple evolutionary routes by which CAM can arise.
Overview of species and experimental set‐up. Left, Tillandsia vanhyningii (Photograph by Michael H.J. Barfuss); right, T. leiboldiana (Photograph by Manuel Lasserus), both at the Botanical Garden of the University of Vienna. Experimental setup, consisting of sampling at four time points that are 6 h apart: 4 and 10 h after dawn (D), plus 4 and 10 h after nightfall (N). Titratable acidity (TA) was quantified for each condition, whereas RNA‐seq data were analysed for the time points indicated with RNA. A separate, additional drought experiment was performed on T. leiboldiana (not shown here) to study net CO2 uptake across diel cycles, to further phenotype its photosynthetic metabolism.
Materials and Methods
2
Experimental Set‐Up and Sampling
2.1
We designed an experiment to test the respective photosynthetic phenotypes of T. vanhyningii and T. leiboldiana in standard (SW) and water‐limited (WL) conditions. The study comprised six clonal accessions per species (Table S1), which were off‐shoots of similar size, as a proxy for age similarity, and were large enough to be considered adult plants able to flower. Accessions from the plant collection of the Botanical Garden of the University of Vienna were placed under identical greenhouse conditions in October 2018. The accessions were allowed to acclimate to the specific greenhouse conditions for 6 weeks, before being placed in a regime of 12 h of light/12 h of dark for 14 days, by complementing daylight with artificial lights between 6 AM and 6 PM. The temperature was maintained at 20°C during nighttime and 22°C during daytime. During these 14 days, three accessions of each species were placed in a water‐limited regime by complete water withholding. The remaining three accessions were kept in standard watering conditions used to maintain the plants in the botanical garden (1× every 2 days). After 14 days, all accessions were sampled every 6 h (four time points: dawn + 4 h [10 AM], dawn + 10 h [4 PM], nightfall + 4 h [10 PM], nightfall + 10 h [4 AM]), resulting in three replicates for the combination of treatments (species × time point × watering regime) (Figure 1). As every accession was sampled four times sequentially, leaf material was collected at each time point by delicately pulling out an entire leaf from the base of the plant to minimise wounding. The tip and base of the leaf were then removed, and the middle part was immediately placed in liquid nitrogen and stored at −80°C.
Titratable Acidity Measurements
2.2
Titratable acidity (TA) was measured following Wai et al. (2019). Frozen leaf material for all four time points, two watering regimes and two species (48 samples in total) was weighed and ground in a precooled mortar to a powder, then added to a 5 mL solution of 80% EtOH. Next, the samples were incubated for 1 h in an 80°C water bath with occasional vortexing. Tubes were centrifuged at 5000 rpm for 5 min, after which the supernatant was transferred and the pellet discarded. The extractions were then let to cool to room temperature. Using a pH526 WTM Sentix mic metre, each sample was measured by titrating 0.1 N NaOH to the solution until reaching a pH of 8.2 exactly and recording the volume of added NaOH. TA in μ‐equivalents per gram of fresh mass (g FM) was calculated using the formula: TA (μ‐equivalents/g FM) = [V NaOH (μL) × C NaOH (μmol/μL)]/m sample (g), where C NaOH is equal to 0.1. Since NaOH is a monovalent base, the number of μ‐equivalents of titratable acid (from both mono‐ and multivalent acids collectively) is equal to the number of micromoles of NaOH consumed during titration.
To test whether TA significantly fluctuated across time within each species and watering regime, we constructed a linear mixed‐effects model (LMM) for each species as follows:
Normality of TA measurements within each species was tested beforehand. In the LMM model, titratable acid represented the dependent variable, while watering regime and time point were fixed effects, and their interaction was included in the model. To control for random variation in TA measurements across replicates of the same time point and watering regime, the model included a random intercept for each replicate. The analysis was performed using the lmer function from R package lme4 (Bates et al. 2015). The goodness of fit of the model was evaluated by plotting a QQ‐Plot and the fitted values to residuals for each species. The significance of fixed effects was evaluated by obtaining a p value using the R package lmerTest (Kuznetsova et al. 2017). Afterwards, a post hoc Tukey HSD test was performed in the R package emmeans to assess the effects of individual watering regimes and time points on TA.
Gas Exchange Measurements of T. leiboldiana
2.3
A separate experiment was performed on T. leiboldiana in October 2025 to assess net CO_2_ uptake progressively throughout a period of water withholding. Leaves of one T. leiboldiana clone were placed in a gas exchange experiment chamber over a period of 8 days where the plant was not watered. The net CO_2_ uptake and other metrics were measured by a WALZ GFS‐3000 system. The ambient temperature during the experiment was maintained between 22°C and 24°C, and the relative humidity at 35% outside the leaf clamp and at 50% inside. Full spectrum LED lights (maximum 340 PPFD) were turned on and off gradually between 6 and 7 AM, and 6 and 7 PM, respectively, resulting in a 12‐h light cycle. At the time of finishing the experiment (9 days without watering), the plant showed clear signs of drought‐stress (yellowing and curling leaves). The plant was then placed under a 14‐day recovery period beginning with full submersion in demineralised water, followed by daily spraying in a greenhouse with same climatic conditions as in the experimental chamber, except for a rH of 70%. Net CO_2_ uptake was measured again at the end of the recovery period. Due to the small size and thick, concave shape of T. vanhyningii leaves, it was not possible to measure gas exchange reliably on this species with this system.
RNA Extraction and Sequencing
2.4
RNA‐seq data was obtained for samples at time points 10 AM (D + 4) and 10 PM (N + 4) for two watering regimes and two species (24 samples in total). Total RNA was extracted using the QIAGEN RNeasy Mini Kit. We evaluated the purity and concentration of the extractions using Nanodrop; RIN and fragmentation profiles were obtained with a BioAnalyzer (Table S1). RNA library preparation and poly‐A capture were performed at the Vienna Biocenter Core Facilities (VBCF) using a NEBNext stranded mRNA kit. The libraries were sequenced on Illumina NovaSeq SP 150 PE.
RNA‐Seq Data Processing
2.5
Raw sequence quality was assessed using FastQC v.0.11.8 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and MultiQC (Ewels et al. 2016). Adapters and low‐quality bases were trimmed off using AdapterRemoval v.2.3.1 (Schubert et al. 2016) with settings –trimns –trimqualities –minquality 20 –trimwindows 12 –minlength 36. The trimmed T. leiboldiana data was aligned to its conspecific reference genome (GenBank accession: GCA_029204045.2, Groot Crego et al. 2024), while the T. vanhyningii data was aligned to the reference genome of its closest relative, T. fasciculata (GenBank accession: GCA_029168755.2; Groot Crego et al. 2024). Plastid phylogenies show that T. fasciculata and the T. ionantha complex, which includes T. vanhyningii (both constitutive, strong CAM), diverged no more than 1.8 Mya (Vera‐Paz et al. 2023). Yardeni et al. (2025) reported 0.008 substitutions per site between the two lineages based on a maximum‐likelihood inference using whole‐genome data. Given this close evolutionary proximity, we consider these species suitable for cross‐species mRNA read alignment. These were performed with STAR v. 2.7.3a (Dobin et al. 2013) in one pass using a GFF file to specify exonic regions. To demonstrate that using a nonconspecific reference genome does not substantially influence our findings, we provide in supplementary a comparison of analyses on T. leiboldiana mapped to its own reference genome, versus to T. fasciculata .
Differential Gene Expression Analysis
2.6
Reads were quantified for paired‐end and reversely stranded reads (options ‐p, ‐s 2) in the subread package FeatureCounts v.2.0.6 (Liao et al. 2014). Counts were generated by counting reads per exon, yet aggregating counts across exons for each gene using FeatureCounts' in‐built meta‐feature function. Differential gene expression (DGE) analyses were performed for each species separately using the Bioconductor package edgeR v.4.2.1 (Robinson et al. 2009).
We built a design matrix incorporating the two investigated treatments (time point and watering regime). Genes where three or fewer samples of the same species had a counts per million (cpm) value smaller than or equal to 1 were filtered out. Library size normalisation was done using edgeR's default trimmed mean of M‐values (TMM) method (Robinson and Oshlack 2010). The Cox‐Reid profile‐adjusted likelihood (CR) method was used to estimate dispersion. Genes were tested for differential expression between time points and between watering regimes separately for each species by fitting a negative binomial general linear model (GLM) with the following contrasts: Day versus Night at SW = Day_SW − Night_SW, Day versus Night at WL = Day_WL − Night_WL, WL vs. SW at day = WL_day − SW_day, WL vs. SW at night = WL_night − SW_night.
To test the fit of these GLMs, we conducted a quasi‐likelihood F‐test. The level of significance was adjusted with the Benjamini–Hochberg False Discovery Rate (FDR) correction to account for multiple testing. Genes were considered differentially expressed (DE) when the FDR‐adjusted p value was below 0.05 and the expression had an absolute log‐fold change of at least 1.5. The number of DE genes for each treatment (time point and watering regime) was then visualised and compared between species using ggplot2.
Overlap in Transcriptomic Response
2.7
To understand the extent of overlap in the transcriptomic response to time and watering regime between species, DE genes of both species were studied on the orthogroup level, using orthology information of T. leiboldiana and T. fasciculata from Groot Crego et al. (2024). Orthogroups containing a DE gene in at least one species are hereafter referred to as DE orthogroups. The percentage of DE orthogroups containing more than one DE gene in a single species was calculated for each comparison to understand the relationship between gene‐level and orthogroup‐level differential expression (reported in Table S8). Venn diagrams were computed to visualise overlaps in DE orthogroups between species using the R package VennDiagram (Chen and Boutros 2011).
We also assessed overlap in DE orthogroups across time‐ and watering‐dependent comparisons, as CAM‐related genes are likely DE in both contexts, particularly in facultative CAM species like T. leiboldiana. The overlap in DE orthogroups between comparisons was visualised with UpSetR (Conway et al. 2017).
In the discussion of results regarding the transcriptomic overlap between species, timepoints and watering regimes, a distinction is made between within‐species comparisons, where the overlap is discussed on the gene level, and between‐species comparisons, where it is discussed on the orthogroup level.
Functional Interpretation of DE Genes
2.8
Gene Ontology (GO) enrichment analyses were performed for both sets of DE genes (between time points and watering regimes) within both species using the Bioconductor package topGO v.2.50.0. (Alexa et al. 2006). Significance of enrichment was evaluated with the Fisher exact test and the algorithm ‘weight01’. The top 100 significant GO terms were then compiled sorted by decreasing significance using a custom‐made script available at the repository https://github.com/cgrootcrego/Tillandsia_CAM‐Drought.
GO enrichment analyses were also performed on each subset of orthogroups representing a section in the Venn diagrams with at least 25 orthogroups. GO enrichments were performed on the DE genes belonging to the orthogroups listed in each Venn diagram intersection. For species‐specific intersections, annotations from the respective species were used (e.g., T. leiboldiana gene annotations for its unique intersections). For overlapping intersections, T. fasciculata annotations were used, assuming overlapping annotations between species within orthogroups.
Complementary to a GO enrichment analysis, we identified genes previously described as related to CAM, starch metabolism, and glycolysis following Groot Crego et al. (2024). Briefly, A. comosus gene IDs related to CAM, starch metabolism or glycolysis compiled by Yardeni et al. (2021) based on a diverse set of studies were used to obtain the corresponding T. fasciculata and T. leiboldiana orthologs using orthology information from Groot Crego et al. (2024). This resulted in a total of 818 and 753 genes, respectively, across 607 orthogroups (Table S2). The orthogroups from this subset that contained DE genes either between time points or watering regimes were then listed by Venn diagram intersection in Tables S8 and S9.
The abundance of RNA‐seq reads for a select group of 53 CAM‐related orthogroups was plotted as a heatmap using the R package ComplexHeatmap (Gu et al. 2016) for all 24 samples. This subset included seven orthogroups containing the pineapple genes putatively described as part of the CAM carbon‐fixation module by Ming et al. (2015), regardless of their differential expression status. The remaining 46 orthogroups comprised aquaporins (5), circadian clock regulators (11), malate transporters (6), photosystem‐related genes (5), vacuolar proton pumps (9), genes related to stomatal movement (4), sugar transporters (5) and one vacuolar invertase (Table S3). Reads were normalised with the TPM (transcripts per million) method using total exonic length of each gene and converted to a log10 scale, then mean‐centred with a z‐score transformation. For homologues of the PEPC and PPCK genes and other core CAM pathway genes, expression in TPM was additionally visualised by line plots.
A separate subset of orthogroups was selected to investigate the drought‐related response to water‐limited conditions in both species, by including all DE genes that were annotated with GO term ‘GO0009414; water deprivation’ or belonged to the drought resistance category as specified by Yardeni et al. (2021). This resulted in 32 orthogroups underlying ABA biosynthesis (1), ABA signalling (4), antioxidation and ROS (reactive oxygen species) scavenging (5), drought tolerance (7), Late Abundance Embryogenesis proteins (8), inositol biosynthesis (3), jasmonic acid (JA) biosynthesis (1), molybdate ion transport (1) and general stress response (2) (Table S4). A heatmap of their expression values in TPM was constructed as detailed above.
Gene Trees of PPC1, PPC2 and PPCK
2.9
We identified homologues of the two subfamilies of PEPC ‐ PPC1 and PPC2 ‐ and of PPCK by running BLASTP on the full set of gene models of A. comosus (V3), Sorghum bicolor (v.3.1.1), Oryza sativa (v.7.0), Zea mays (RefGen v.4), Asparagus officinalis (v.1.1) and Arabidopsis thaliana (TAIR10) available on Phytozome (https://phytozome‐next.jgi.doe.gov/), with the PEPC and PPCK gene models of T. fasciculata and T. leiboldiana as query (Groot Crego et al. 2024). Hits were filtered based on identity and target length for each respective gene. For PPC1, we selected all targets longer than 500 aa and with more than 80% identity, resulting in 56 protein sequences across the seven species. For PPC2, the identity threshold was lowered to 70%, resulting in eight sequences. For PPCK, we selected all sequences with an identity greater than 60% regardless of target length, resulting in 17 sequences.
Protein sequences were aligned using mafft v.7.520 (Katoh et al. 2002) with options –op 1.5 –ep 0.1 –maxiterate 1000 –localpair. Gene trees were constructed using IQTREE v.2.2.5 (Minh et al. 2020) with options ‐B 1000 ‐T AUTO. The resulting consensus tree was used for visualisation, collapsing branches leading to Poaceae genes. Nomenclature of PEPC homologues follows Deng et al. (2016) and Heyduk et al. (2022).
Results
3
Titratable Acidity Measurements
3.1
In standard conditions (SW), TA displayed a CAM‐like diel curve in T. vanhyningii, with overnight acid accumulation followed by daytime depletion (Figure 2a, Table S5). Our LMM model showed that TA significantly fluctuated between time points both under standard and water‐limited (WL) conditions, yet there was no significant difference in TA values between watering regimes at the same time point (Figure 2, Figure S1, Table S6). In T. leiboldiana, TA did not significantly accumulate overnight under SW (Figure 2b). Under WL, TA values between D + 10 and N + 10 were significantly different, indicating a slight fluctuation in leaf acidity (Figure 2b). T. leiboldiana exhibited a markedly different temporal curve in TA accumulation between SW and WL, as TA values between watering regimes at each time point was significantly different in all but one time point (Figure 2a, Table S6). This indicates a weak, facultative CAM cycle in T. leiboldiana, whereas T. vanhyningii exhibited temporal TA signatures that are characteristic of a strong, constitutive CAM plant.
*Titratable acidity (TA) measurements of Tillandsia vanhyningii and T. leiboldiana. (a) TA measured at four time points across a diel cycle under two watering regimes. The significance level of within‐species, between‐watering regime comparisons at a given time point from an LMM model are shown with following thresholds: ns—p value > 0.05; *p value < 0.05; **p value < 0.01; **p value < 1E‐03. (b) Difference in adjusted means of TA values between each pair of timepoints, separated for each species. Each dot represents the difference between adjusted means of TA values between two timepoints within a species and watering regime, as estimated in the post hoc Tukey HSD test. The first four timepoint comparisons are between adjacent timepoints that are 6 h apart in a 24‐h cycle, whereas the last two represent differences in mean TA between timepoints that were 12 h apart. The shape of each dot indicates the significance level of the difference in mean TA and is shown with the same thresholds as in (a). D, dawn; FM, fresh mass; N, nightfall; SW, standard watering conditions; WL, water‐limited conditions.
Gas Exchange in T. leiboldiana
3.2
At the start of the water‐withholding experiment, T. leiboldiana showed a C_3_‐like pattern of CO_2_ assimilation occurring principally during the day (Figure 3a, Table S13). However, throughout Day 2 and 3, net daytime CO_2_ uptake decreased in a step‐wise manner till null through most of Day 3. In Days 1–3, there was no night‐time net CO_2_ uptake. On Day 4, a metabolic shift occurred towards a weak CAM‐like pattern of night‐time CO_2_ assimilation. From this point onwards, night‐time net CO_2_ uptake became visible and increased gradually as the experiment progressed, with the highest values on Day 8. T. leiboldiana also exhibited net CO_2_ uptake during the day, but this was consistently higher at night. After a 14‐day recovery, the plant showed day‐time net CO_2_ uptake and none during the night, an indication that it had fully shifted back to a C_3_‐like photosynthetic metabolism (Figure 3b, Table S14).
Gas exchange measurements of Tillandsia leiboldiana under a progression from well‐watered to water‐limited conditions, followed by recovery. Net CO2 uptake (blue line, left Y‐axis) and photosynthetically active radiation (PAR; red line, right Y‐axis) in leaves of T. leiboldiana. For visualisation purposes, exponential smoothing was applied to the net assimilation curves. When the value of net CO2 uptake is positive (above horizontal line), the plant is assimilating more CO2 than it is leaking CO2, the opposite is true for negative values. A 24‐h cycle is divided into a light period (white background) and a dark period (grey background). (a) Measurements throughout an 8‐day period without watering. (b) Measurements across one 24‐h cycle after a 14‐day recovery.
Differential Gene Expression
3.3
Unique mapping rates were high for both species, but T. leiboldiana's reads mapped better to its own genome than T. vanhyningii's reads to the T. fasciculata genome (Figure S2). On average, T. leiboldiana had 16.3% more uniquely assigned reads per sample than T. vanhyningii (Table S7). However, T. vanhyningii had 8.7% more genes with non‐zero counts compared to T. leiboldiana. A comparison of downstream results when mapping T. leiboldiana to its own reference genome, versus to T. fasciculata , shows minimal effect (see Supporting Information). T. leiboldiana is approximately twice as distantly related to T. fasciculata as T. vanhyningii is. Therefore, we expect a much smaller effect using T. vanhyningii alignments to T. fasciculata on the results presented in this study.
In T. vanhyningii, variance in read counts was largely explained by time point, while watering regime was partially resolved in the second component of the PCA (Figure S3). In T. leiboldiana, the first PC separated samples by watering regime, and the second component showed a clear split by time point. Overall, T. leiboldiana showed clearer transcriptomic distinctions by time point and watering regime than T. vanhyningii (Figure S3).
DGE analyses revealed widely different counts of DE genes along the time and watering regime axes. In both species, many genes were DE between day and night, regardless of watering regime (Figure 4). Under SW, the number of upregulated genes in T. vanhyningii was comparable between day and night, with only six more genes up‐regulated at night than at day, but T. leiboldiana had 2.4× more genes upregulated at day than at night. Under WL, more genes were upregulated at night in both species compared to SW. While this pattern is noticeable in T. vanhyningii, with 57% of DE genes upregulated at night under WL compared to 50% under SW, it is especially pronounced in T. leiboldiana (46% of genes upregulated at night under WL compared to 30% under SW).
Counts of up‐ and downregulated DE genes (FDR‐adjusted p value < 0.05, |Log(FC)| ≥ 1.5) in four comparisons across different background conditions. The first four bars show counts of time‐dependent DE genes in two watering regimes, while the last four bars show counts of genes that are DE between watering regimes at two time points. Left‐hand symbols indicate the background condition, for example, the first two bars are counts of time‐dependent DE genes under standard conditions. Left‐hand bars show counts of genes that are upregulated at day compared to night for time‐dependent comparisons, and in water‐limited compared to standard conditions for comparisons of watering regimes. Right‐hand bars show counts of genes that are upregulated at night compared to daytime for time‐dependent comparisons, and in standard watering compared to water‐limited conditions for comparisons of watering regimes. The DE analyses were conducted in each species separately. A dark‐grey background indicates night‐time, while lighter bar colours indicate water‐limited conditions.
In comparison with time‐dependent analyses, fewer genes were DE between watering regimes, except in T. leiboldiana at night. Contrary to T. leiboldiana, T. vanhyningii showed a very low number of DE genes between watering regimes, independent of time point, suggesting that most temporal expression patterns detected in the time‐dependent DGE analysis were constitutive.
Compared to T. vanhyningii, T. leiboldiana manifested a more noticeable response to water withholding, especially at night. During the day, most genes appeared downregulated under WL compared to SW (Figure 4). The total number of DE genes between watering conditions at night was 2.9× larger than at day, which is a contrary trend to T. vanhyningii, where 2.5× more genes were DE between watering regimes at day than at night.
Overlap in Transcriptomic Response
3.4
The time‐dependent transcriptomic response showed large heterogeneity across species and watering regimes, with 50% of all orthogroups with one or more DE genes (hereafter DE orthogroups) being unique to one species and watering regime (Figure 5a, sections a, b, c and d, Table S8). Only 5.9% of all time‐dependent DE orthogroups were shared between all groups (section abcd). These were enriched for photosynthesis and circadian clock‐related functions (Table S9).
Overlap of DE orthogroups between species, time points and watering regimes. (a) Venn diagram showing the number of shared time‐dependent DE orthogroups (i.e., day versus night) per combined group of species and watering condition. Example: Section a shows the number of orthogroups that contained genes that were exclusively DE between day and night under standard watering conditions in Tillandsia vanhyningii, while section ad shows the number of orthogroups with a gene that is DE both in T. leiboldiana and in T. vanhyningii under standard watering conditions. (b) Venn diagram showing the number of shared DE orthogroups between watering regimes (i.e., water‐limited versus standard conditions) per combined group of species and time point. Letters inside sections show the overlap between the four groups. Example: Section a shows the number of orthogroups that contained genes that were exclusively DE between watering regimes at day‐time in T. vanhyningii, while section ad shows the number of orthogroups with a gene that is DE both in T. leiboldiana and in T. vanhyningii at day‐time.
The largest overlap was found within species, between conditions, as only 19% of all time‐dependent DE genes were shared between species (Table S8, intersections of a and b, and c and d). In T. leiboldiana, genes that were uniquely time‐dependent DE under WL and not under SW (Figure 5a, Table S9, sections d, bd, abd and ad) were enriched for responses to water deprivation, abiotic stress and abscisic acid (ABA), while genes that were only time‐dependent DE under SW (Figure 5a, Table S9, sections c, bc, abc and ac) lacked such enrichments, indicating that T. leiboldiana experienced drought‐induced stress under WL (Table S9). Conversely, functions related to water deprivation, response to abiotic stress or ABA were enriched in time‐dependent DE genes shared between watering regimes in T. vanhyningii (Figure 5a, Table S9, sections ab, abc, abcd and abd).
Under WL, more time‐dependent DE orthogroups were shared between species than under SW (17% vs. 14%, Table S8), with enriched functions related to stomatal movement, ABA, and starch metabolism (Figure 5a, Table S9, sections ad, acd and abd).
The response to watering regime was largely nonshared between groups, with 75% of DE orthogroups being unique to one species and time point (Figure 5b, sections a, b, c, and d) and only two DE orthogroups being shared among all groups (Figure 5b, section abcd). T. vanhyningii showed little overlap between time points in the response to WL (8%), while in T. leiboldiana 24% of genes are DE between watering regimes at both time points (Table S8). Due to the small number of DE genes in T. vanhyningii between watering regimes, the overlap between species was very small (2.8%).
In summary, both the overall transcriptomic response to day‐night changes and to watering regimes showed little overlap between species, yet diel expression between species seemed to overlap slightly more under water‐limited than under standard conditions, indicating an increase in the resemblance of the two species' metabolism under water‐limited conditions.
CAM‐Related Transcriptomic Overlap and Differences Between Species, Time Points and Conditions
3.5
Overlap Between Species
3.5.1
In T. leiboldiana, CAM‐related gene expression is likely to be associated both with differential expression between time points and between watering regimes, as it activated a CAM cycle under water‐limited conditions (Figures 2 and 3). When considering all DE orthogroups summed up across species, time points and watering regimes, 516 contained genes that were both DE between time points and watering regimes across species.
Of these 516 orthogroups, the largest group was exclusively DE in T. leiboldiana, both between time points and watering regimes (Figure S4), and included orthogroups encoding for aquaporin NIP5‐1, an aluminium‐activated malate transporter (ALMT) 12 and a sugar transporter 5. The second largest group consisted of 132 orthogroups (26%) which displayed time‐dependent expression in T. vanhyningii while being DE between watering regimes in T. leiboldiana, and included 12 previously described CAM‐related orthogroups. Most notable were orthogroups encoding for the sugar and malate transporters SWEET14 and ALMT 9. SWEET14 showed elevated expression under WL in both species, yet it seemed upregulated at night in T. leiboldiana, while in T. vanhyningii its highest expression was during the day (Figure 6). ALMT 9 was more highly expressed in T. leiboldiana under SW, and showed lower expression under WL. T. vanhyningii's expression of ALMT 9's ortholog was comparatively lower to T. leiboldiana's across all time points and treatments (Figure 6).
Transcript abundance in ‘transcripts per million’ (TPM) for a subset of single‐copy CAM‐related orthogroups (listed in Table S3), in Tillandsia vanhyningii (left panel, using the T. fasciculata ortholog) and T. leiboldiana (right panel, using the T. leiboldiana ortholog) at two time points and watering regimes. TPM values are transformed to a logarithmic scale and centred as a z‐score. Displayed genes are organised in functional categories and gene name colour reflects their differential expression status: Blue = DE between watering regimes, green = DE between time points, black is DE in both comparisons and grey = not DE. The letters in superscript indicate the intersection of the Venn diagram that each orthogroup belongs to. Schematics of Venn diagrams in the bottom left indicate what section each letter refers to. When an orthogroup is DE in both comparisons, the first superscripted letter indicates the intersection of the time‐dependent Venn diagram (Figure 5a), and the second superscripted letter, separated by a comma, indicates the intersection of the Venn diagram for watering regimes (Figure 5b). ALMT, aluminium‐dependent malate transporter; ARF, auxin response factor; bZIP, basic region leucine zipper; LEA, late embryogenesis abundant; LHCB, chlorophyll a–b binding protein; LHY, late elongated hypocotyl; NAD‐MDH, NAD‐dependent malate dehydrogenase; NAD‐ME, NAD‐dependent malic enzyme; NIP, nodulin 26‐like intrinsic protein; OEE, oxygen evolving enhancer; PEPCK, phosphoenolpyruvate carboxykinase; PIF, phytochrome interacting factor; PIP, plasma membrane intrinsic protein; PPC, phosphoenolpyruvate carboxylase; PPCK, phosphoenolpyruvate carboxylase kinase; PPDK, pyruvate‐phosphate dikinase; PPDK‐reg. prot, PPDK‐regulatory protein; RUP, repressor of UV‐B photomorphogenesis; SLAC1, slow‐type anion channel 1; STP, sugar transporter protein; SUT1, sucrose transporter 1; TDT, tonoplast dicarboxylate transporter; TIP, tonoplast intrinsic proteins; α‐CA, alpha carbonic anhydrase; β‐CA, beta carbonic anhydrase.
One hundred twenty‐seven orthogroups (25%) showed a time‐dependent expression in both species but were only DE between watering conditions in T. leiboldiana (Figure S4). Among these were 18 CAM‐related orthogroups, yet only seven contained genes that were up‐regulated under water‐limited conditions in T. leiboldiana. This included an orthogroup encoding for PEPC kinase (PPCK2), which is crucial for CAM photosynthesis, as it activates the main carboxylating enzyme phosphoenolpyruvate carboxylase (PEPC, see Section 3.5.4 for detailed discussion). Additionally, an orthogroup containing the drought resistance gene YLS9 (Müller et al. 2017), and three orthogroups encoding for enzymes associated with glycolysis and gluconeogenesis were part of this group. YLS9 showed time‐dependent expression in T. vanhyningii in both watering regimes, while in T. leiboldiana, it acquired weakly time‐dependent expression under WL (Figure 6). Together with PPCK2, YLS9 was one of the few CAM‐related genes that exhibited an expression profile in T. leiboldiana that more closely resembled that of T. vanhyningii under WL.
We also searched for CAM‐related orthogroups among the larger set of DE orthogroups called across time and watering regime separately. In T. vanhyningii and T. leiboldiana, 97 and 75 time‐dependent DE CAM‐related orthogroups were identified, respectively, revealing functions like circadian clock, glycolysis, stomatal movement, drought resistance and water/sugar transport (Tables S5 and S8). Overlap between species and watering regimes of these orthogroups followed a similar trend as the full set (51% of genes unique to a species and watering regime). Twelve CAM‐related DE orthogroups shared across all groups were linked to the circadian clock, photosynthesis and photosystem activity (Table S10, section abcd). T. vanhyningii shared 46% of CAM DE genes between watering conditions, while T. leiboldiana shared 36%, indicating a stronger metabolic shift under WL. The proportion of time‐dependent CAM DE orthogroups shared between species was slightly higher under WL (25%) than under SW (22%), yet remains a minority of the total set of CAM‐related DE genes.
For orthogroups that exhibited DE between watering regimes, six and 75 CAM‐related orthogroups were present in T. vanhyningii and T. leiboldiana, respectively (Tables S8 and S11). Only one orthogroup overlapped between species, whose orthologs encode for a probable linoleate 9S‐lipoxygenase enzyme found differentially expressed between C_3_ and CAM Tillandsia species in a previous study (De La Harpe et al. 2020). Lipoxygenases are involved in the biosynthesis of jasmonates, which play an important role in response to various biotic and abiotic stressors (Bachmann et al. 2002).
CAM‐Related Expression in Tillandsia vanhyningii
3.5.2
In SW, T. vanhyningii showed time‐dependent differential expression of core CAM gene PPCK and of genes underlying stomatal movement, vacuolar proton pumps, malate transport, sugar transport and glycolysis (Table S10). Several genes encoding for Late Embryogenesis Abundant (LEA) proteins linked to drought tolerance (Arumingtyas et al. 2013; Magwanga et al. 2018) also showed diel expression in SW.
Forty‐six CAM‐related time‐dependent DE genes overlapped between SW and WL (Table S8), including genes encoding for PPCK, circadian clock regulator REVEILLE1, a beta‐amylase (starch breakdown), two auxin response factors (ARF, linked to stomatal function; Bouzroud et al. 2020), two aquaporins and multiple LEA proteins (Table S10), indicating that genes from a wide set of CAM‐related functions were expressed similarly between conditions. AVP1, a secondary proton pump that aids in the vacuolar transport of malate, was only DE in SW. A gene encoding for a chloroplastic malate dehydrogenase (MDH) showed diel expression in both conditions, though it was only significantly time‐dependent under WL (Figure 6). It was up‐regulated during the day, suggesting a role in the day‐time decarboxylation module of CAM. Its corresponding ortholog in T. leiboldiana showed no time‐dependent expression in any condition. We also observed a shift in time‐dependent expression between conditions among genes encoding for sugar transporters. While sugar transporter SWEET16 showed significant time‐dependent expression exclusively under SW, SWEET14 was only DE under WL (Figure 6). Both genes were up‐regulated during the day when displaying time‐dependent expression.
Very few CAM‐related genes were DE in T. vanhyningii between watering regimes. They included several genes encoding enzymes associated with glycolysis and gluconeogenesis and genes previously described as DE between C_3_ and CAM Tillandsia (De La Harpe et al. 2020). None of these genes were part of the core CAM pathway, and while genes encoding many core CAM enzymes did show a slight increase in expression under water‐limited conditions (β‐carbonic anhydrase, NAD‐dependent MDH, NAD‐ME, Figure 6), their expression was similar overall between watering regimes.
CAM‐Related Expression in Tillandsia leiboldiana
3.5.3
Under SW, time‐dependent CAM‐related DE genes were associated with functions that are typically expected to show temporal oscillation in plants that do not perform CAM (glycolysis and gluconeogenesis, sugar transport, circadian clock regulation). Under WL, 29 DE genes overlapped with standard conditions (Table S10), yet many additional CAM‐related genes became time‐dependent DE, involved in malate transport, stomatal movement, sugar metabolism, water transport and starch metabolism. Genes encoding Alpha‐carbonic anhydrase 1 (α‐CA) and phosphoenolpyruvate carboxylase kinase (PPCK) were upregulated at night. SLAC1 and two auxin response factors related to stomatal function were also time‐dependent DE under WL (Table S10, Figure 6). These findings further emphasise that a CAM cycle was activated in T. leiboldiana under WL.
Furthermore, most CAM‐related genes that were DE between watering regimes could be found in T. leiboldiana at night (Venn diagram section c) and included genes encoding the core CAM enzymes PEPC and PPCK (Table S11). Several genes encoding malate transporters were also in this group, though they were generally downregulated under WL, except for ALMT10 (Figure 6). Sugar transport genes were upregulated under WL, such as genes encoding sugar transporter proteins (STPs) and Sucrose Transporter 1 (SUT1), whose ortholog also showed elevated expression under WL in T. vanhyningii. Genes related to stomatal function were generally more highly expressed under WL but displayed varying behaviour in terms of the time of expression (Figure 6). While ARF17 and SLAC1 were upregulated during the day, ARF19 and TRAB1 were upregulated during the night.
Diel Expression of PPC and PPCK Homologues, and Other CAM‐Related Gene Families, Between Species and Watering Regimes
3.5.4
The transcript abundances of PEPC and PPCK homologues revealed a variety of expression patterns. Gene Tfasc_v1.16311‐RA, which is the ortholog of PPC1‐M1 (Figure S5a) and was previously identified as the CAM‐related fluctuating copy of PEPC in T. fasciculata (Groot Crego et al. 2024) did not show clear diel expression changes in T. vanhyningii (Figure 7). In T. leiboldiana a slight increase was noticeable in its ortholog at night but was not significantly DE in either condition. However, PPC1‐M1 was noticeably more highly expressed in T. vanhyningii compared to T. leiboldiana. The remaining known PEPC genes were generally less expressed than PPC1‐M1 in both species but appeared to become activated under WL, namely in T. vanhyningii for PPC1‐M2 (not significant), and in T. leiboldiana for PPC2.
Transcript abundance in ‘transcripts per million’ (TPM) for all homologous genes encoding for PEPC and PEPC kinase (PPCK) in Tillandsia vanhyningii (left panel) and T. leiboldiana (right panel) at two time points and watering regimes. TPM values are transformed to a logarithmic scale. The colour of gene names reflects their differential expression status: Blue = DE between watering regimes, green = DE between time points, black is DE in both comparisons and grey = not DE. The superscripted letters after a gene name indicate the intersection of the Venn diagram this gene belonged to.
For annotated PEPC kinases (PPCKs), a diel expression was strongest in the PPCK2 homologue (Figure S6) in T. vanhyningii under both watering regimes. While its ortholog was not strongly expressed at night under SW in T. leiboldiana, it became highly expressed under WL (Figure 7). PPCK1 showed mild diel fluctuation in T. vanhyningii, but its ortholog was overall lowly expressed in T. leiboldiana.
The NAD‐dependent malate dehydrogenase (NAD‐MDH, OG0001548), previously identified as CAM‐related (Groot Crego et al. 2024), is crucial for night‐time carboxylation in CAM. This gene has two copies in A. comosus , T. fasciculata and T. leiboldiana (Figure S7a). Neither copy showed significant differential expression in this study. One copy had increased night expression, while the other was higher during the day (Figure S7b), with both showing slightly elevated expression under WL in both species. Despite similar time‐dependent expression profiles, the constitutive CAM plant exhibited substantially higher expression of this and other genes encoding for core CAM enzymes, including V‐ATPase, PEPCK, NAD‐ME, PPDK and PPDK‐regulatory protein (Figure 6, Figures S8–S10).
Description of the General Response to Water Withholding
3.6
GO enrichment of watering condition‐dependent DE genes in T. leiboldiana revealed a pronounced drought response, with the top upregulated genes linked to water deprivation and ABA signalling. Processes that were downregulated under WL pointed at a general shutdown of metabolism, such as photosynthesis, transport and biosynthetic processes (Figure 8a, Table S12). In contrast, T. vanhyningii showed minimal drought‐related responses, with no significant changes in regulation of photosynthesis or water deprivation genes. Yet, a few enriched GO terms suggested a limited drought response, including inositol, jasmonate and molybdate metabolism (Figure 8b, Table S12).
Summary of significantly enriched GO terms representing a biological process (BP) in DE genes between watering regimes in Tillandsia leiboldiana (a) and T. vanhyningii (b). GO term enrichment at day and night were combined for each bar plot. The bars show the significance of each enriched GO term as the p value adjusted for multiple testing, on a logarithmic scale. The colour gradient indicates up‐ and downregulation of the genes underlying each GO term, as a z‐score of the fold change on a logarithmic scale. The GO terms are shown in descending order from most upregulated in water‐limited conditions to most upregulated in standard conditions.
We detected 32 orthogroups related to drought response and water deprivation that were differentially expressed across species, watering regimes and time points. Eight of these were upregulated under WL in T. leiboldiana compared to SW, and four were downregulated (Figure S12). The upregulated orthogroups encoded for three LEA proteins, a phosphatase 2C, which is involved in ABA signalling (Jung et al. 2020), a galactinol‐sucrose galactosyltransferase, which is associated with ROS scavenging (Nishizawa et al. 2008), an annexin which also confers oxidative protection under stressful conditions (Konopka‐Postupolska et al. 2009), and the protein YLS9 (see Section 3.5.1).
In T. vanhyningii, only five orthogroups related to drought‐induced stress were differentially expressed between watering regimes, of which three were upregulated under WL: an inositol‐3‐phosphate synthase, a 12‐oxophytodienoate reductase involved in jasmonic acid biosynthesis and a molybdate ion transporter. In contrast to the low number of drought‐related genes that were DE between watering regimes in T. vanhyningii, 18 were DE between time points. However, out of these, 16 already exhibited time‐dependent expression under SW. Several orthologs also showed overall higher expression in T. vanhyningii than in T. leiboldiana. Genes encoding LEA proteins 5 and 8, for example, were more highly expressed in T. vanhyningii and showed diel oscillation in both watering regimes, with their expression increasing at night (Figure S12). Orthologous genes related to ABA signalling exhibited time‐dependent expression in T. vanhyningii under both watering regimes (except for phosphatase 2C), and in some cases also in T. leiboldiana (Figure S12). A 9‐cis‐epoxycarotenoid dioxygenase (NCED), involved in ABA biosynthesis (Wan et al. 2011), also showed diel expression changes in T. vanhyningii under SW. These findings suggest that high and/or diel expression of drought‐related genes may be an inherent feature of constitutive CAM plants.
Discussion
4
CAM is associated with diversification in many plant families and is considered a key innovation in Bromeliaceae (Benzing 2000; Ogburn and Edwards 2009; Quezada and Gianoli 2011; Silvera et al. 2009). However, its evolutionary trajectory and the degree of genetic overlap between distinct CAM forms are not well understood (Edwards 2023). The subgenus Tillandsia displays a wide spectrum of CAM phenotypes (Crayn et al. 2015; De La Harpe et al. 2020), but this diversity has so far been characterised only partially, mostly through carbon isotope ratios and transcriptomics under standard conditions. By comparing physiological and temporal transcriptomic responses to water limitation in two Tillandsia species, we aimed to explore the phenotypic range of CAM in this group and assess to what extent facultative and constitutive CAM rely on shared transcriptomic mechanisms.
Tillandsia leiboldiana Is a Facultative CAM Plant
4.1
Our TA, gas exchange and DGE analyses point at the presence of a latent CAM cycle in T. leiboldiana that is activated under drought‐induced stress. First, water‐withholding induced weak, yet significant night‐time TA accumulation, absent in well‐watered plants. Nocturnal acidification is a diagnostic feature of CAM (Winter and Smith 2022) and has previously been used to identify facultative CAM in Clusia (Borland et al. 1992), Mesembryanthemum (Winter and Ziegler 1992) and Agave (Heyduk et al. 2016).
Second, we observed a switch to net nocturnal CO_2_ uptake after 4 days of withholding water, followed by full recovery of day‐time uptake after re‐watering (Figure 3b). Nocturnal net CO_2_ uptake is another hallmark of CAM, indicating open stomata and storage of CO_2_ in the vacuole in the form of malate (Lüttge 1987), which aligns with our TA measurements. While the amplitude of nocturnal CO_2_ uptake was low, this pattern is comparable to other facultative CAM plants under drought. The CO_2_ uptake profiles in drought‐stressed Clusia pratensis, Calandrinia polyandra, Talinum triangulare and Cistanthe cachinalensis all show modest night‐time rates, and some daytime net assimilation remains (Winter and Holtum 2014; Chomentowska et al. 2025).
Transcriptomic signatures further support altered stomata behaviour for T. leiboldiana. Differential expression of stomatal function genes between watering regimes suggests altered stomatal behaviour under WL, though this may reflect a general drought‐stress response rather than CAM‐specific activity. For instance, TRAB1, involved in ABA signalling during drought stress (Hobo et al. 1999; Zhu 2002), is upregulated in T. leiboldiana under WL at night (Figure 6). However, the guard cell anion channel SLAC1 controlling stomatal closure (Vahisalu et al. 2008) is downregulated at night under WL (Figure 6), potentially indicating that stomata open at night.
Third, night‐time upregulation of PPCK2 in T. leiboldiana under WL indicates induction of the phosphorylation‐based regulatory machinery necessary for nocturnal PEPC activity—a pattern that is entirely absent in SW.
Carbon isotope ratios reported for T. leiboldiana (δ^13^C = −28.0 to −31.3; Crayn et al. 2015; De La Harpe et al. 2020) would typically lead to its classification as a C_3_ plant. However, isotope values can mask facultative CAM activity, as seen in other bromeliads such as Guzmania and Ronnbergia (Pierce et al. 2002). Our integrative study therefore reveals that T. leiboldiana is not strictly C_3_, but a weak, facultative CAM species that performs C_3_ photosynthesis under benign conditions while acquiring weak CAM function under drought. Minimal expression changes in core CAM pathway genes between watering regimes further suggest that high transcript abundance of these genes may precondition facultative CAM expression in Tillandsia.
This raises an evolutionary question: What is the point of origin of CAM in the subgenus Tillandsia? Ancestral state reconstructions using carbon‐isotope data support a C_3_ ancestor for the subgenus (De La Harpe et al. 2020) and multiple independent origins of CAM within the clade. However, the presence of inducible CAM in T. leiboldiana opens an alternative possibility—that at least parts of the CAM cycle may have been present prior to the divergence of the lineages leading to T. leiboldiana and T. vanhyningii. Although our pairwise design cannot rule out independent origins of facultative and constitutive CAM, it highlights the need for broader comparative datasets.
T. vanhyningii Shows Minimal Physiological and Transcriptomic Response to Water Withholding
4.2
Our physiological and transcriptomic analyses indicate that the quantitative response to water‐withholding was far greater in T. leiboldiana than in T. vanhyningii (Figures 2, 3, 4), suggesting that T. leiboldiana's metabolism was affected by water deprivation, while T. vanhyningii's seemed largely the same between conditions and therefore not substantially stressed. TA accumulated in a similar manner between watering regimes in T. vanhyningii, though accumulation occurred faster overnight under water‐limited conditions. This aligns with studies on T. ionantha, which showed minimal physiological changes after 6 months of water withholding, attributed to daytime stomatal closure, low stomatal density and efficient atmospheric moisture absorption via its dense trichome layer (Nowak and Martin 1997; Ohrui et al. 2007).
A very small number of genes were found to be DE between watering regimes in T. vanhyningii, while a prominent response was detected in T. leiboldiana. Upregulated enriched GO terms under WL in T. vanhyningii included metabolism of inositol, jasmonate, and molybdate (Figure 8b, Table S12), all known to play a role in drought‐stress response (Li et al. 2022; Loewus and Murthy 2000; Wasternack 2014). However, the expression of many genes related to drought and abiotic stress response was time‐dependent or more highly expressed in T. vanhyningii compared to T. leiboldiana under SW, including ABA biosynthesis and signalling genes or LEA protein genes. A link between constitutive CAM phenotypes and temporal expression of ABA‐related genes was observed by De La Harpe et al. (2020). ABA signalling plays an important role in stomatal closure as a response to water deprivation (Assmann and Jegla 2016; Christmann et al. 2006; Wilkinson and Davies 2002), but has also been associated with regulation of CAM photosynthesis (Mioto and Mercier 2013; Ting 1981). Our findings corroborate that time‐dependent expression of ABA‐related genes is characteristic of constitutive CAM phenotypes and does not necessarily indicate that T. vanhyningii was drought‐stressed in our experiment.
Regarding the CAM‐related response to WL, our findings showed that the expression of CAM‐related genes in T. vanhyningii remained largely similar. In some cases, however, such as sugar transporters SWEET16 and SWEET14, different genes performing the same functions were recruited for CAM under varying watering regimes.
Atmospheric or ‘grey’ Tillandsia such as T. vanhyningii represent the most extreme form of adaptation to an epiphytic lifestyle among Bromeliaceae, having largely lost the ability to absorb water with roots (Benzing 2000; Crayn et al. 2004). They occur in arid habitats and are all constitutive CAM plants, suggesting that CAM is a key innovation in this group (Crayn et al. 2004). Our findings confirm that atmospheric Tillandsia are drought‐tolerant, as previously shown physiologically (Castillo et al. 2016; Martin 1994; Nowak et al. 1997; Stiles and Martin 1996), but now supported at the transcriptomic level, as T. vanhyningii was only mildly affected by water‐limited conditions. However, while this could be attributable to T. vanhyningii's constitutive CAM phenotype, other typical characteristics of atmospheric Tillandsia such as a dense trichome layer, need to also be considered as important contributors to their drought‐tolerant nature.
Shared Versus Non‐Shared Transcriptional Architecture of Photosynthetic Metabolism
4.3
Despite the relatively recent divergence of T. vanhyningii and T. leiboldiana (Barfuss et al. 2016; Mendoza et al. 2017; Yardeni et al. 2025), the two species show clear physiological and transcriptomic differentiation consistent with distinct photosynthetic phenotypes. Only limited overlap in time‐dependent and watering‐dependent DE genes was found between species, mirroring observations in other CAM lineages such as Agavoideae (Heyduk et al. 2022).
Most CAM‐related, time‐dependent DE genes showed diel expression cycling in only one species (Figure 9). Although T. leiboldiana seemingly activated CAM under WL, the underlying transcriptomic architecture differed substantially from the constitutive CAM pattern of T. vanhyningii. Notable exceptions were PPCK2 (Figures 7 and 9) and the drought‐associated protein YLS9, which showed a similar time‐dependent expression both in T. vanhyningii and drought‐stressed T. leiboldiana. This restricted overlap suggests that the facultative and constitutive CAM cycles of these species rely on largely distinct transcriptional solutions, rather than expression‐level modulation of a common pathway.
Schematic of the CAM pathway (adapted from Groot Crego et al. 2024), with an overview of CAM‐related gene expression in standard watering and water‐limiting conditions in Tillandsia vanhyningii and T. leiboldiana. Enzymes are shown in boxes, while pathway products are bold outside boxes. Enzymatic members of CAM metabolism pathways are shown in orange, while stomatal and circadian regulators are highlighted in yellow. Stomatal regulators are on the left outside the cell and circadian regulators on the right. Symbols indicate whether genes were DE in this study in timewise comparisons (clock symbol) or between watering regimes (watering can symbol). The colour of the symbol indicates in which species the differential expression was detected, while a gradient of both colours indicates DE in both. ‘Different homologues’ indicates that the species show CAM‐like expression in distinct homologues of a gene family underlying the enzyme. ‘Shared expression’ indicates high similarity in expression pattern between species across conditions. ARF, auxin response factor; CA, carbonic anhydrase; GI, protein GIGANTEA; LHY, late elongated hypocotyl; MDH, malate dehydrogenase; PEPC, phosphoenolpyruvate carboxylase; PEPCK, phosphoenolpyruvate carboxykinase; PPCK, phosphoenolpyruvate carboxylase kinase; RVE1, REVEILLE 1; SLAC1, slow‐type anion channel 1; SUT1, sucrose transporter 1.
This is supported by the distinct expression patterns of gene family members of PEPC, a core CAM gene. Under WL, the homologous copy of PPC2 appeared upregulated in T. leiboldiana, while being lowly expressed under all conditions in T. vanhyningii. Because none of the PEPC homologues are time‐dependently expressed in this experiment, likely due to the time of sampling not being synchronous with the peak expression times of these genes, it is difficult to confirm that PPC2 is active in CAM in T. leiboldiana under drought‐stress. Both in T. fasciculata and in Agave, PEPC gene expression peaks at dusk, while samples were taken 4 h after dusk in this study (Abraham et al. 2016; Groot Crego et al. 2024). While in the majority of CAM plant lineages, including T. fasciculata (Groot Crego et al. 2024) and A. comosus (Ming et al. 2015), the PEPC homologue recruited into the CAM pathway is PPC1, recruitment of PPC2 has been suggested in Yucca aloifolia (Heyduk et al. 2022) and in Isoetes taiwanensis (Wickell et al. 2021). This supports the view that PEPC paralog recruitment can evolve flexibly within CAM lineages.
Another core CAM gene that showed increased time‐dependent expression in T. leiboldiana under WL was α‐carbonic anhydrase (α‐CA, Figure S11). Carbonic anhydrases convert CO_2_ to HCO_3_ ^−^, the first step of carbon assimilation in CAM photosynthesis. In pineapple, however, β‐CA, rather than α‐CA, is modulated in CAM (Ming et al. 2015), and T. vanhyningii does not share the time‐dependent expression of this gene. In several Yucca species, on the other hand, α‐CA did show time‐dependent expression under water‐limited conditions (Heyduk et al. 2019). The CAM‐like expression of PPC2 and α‐CA in T. leiboldiana suggests species‐specific recruitment of different genes and homologues into the CAM pathway (Figure 9).
Except for PPCK, PPC2 and α‐CA, the core (de)carboxylation module seems relatively conserved between species. However, other CAM‐associated processes (Figure 9), such as malate transport, sugar transport, stomatal regulation, water channels and starch metabolism, exhibited substantial transcriptomic divergence between species. Overall, the realised CAM phenotypes are underpinned by a complex mixture of conserved patterns at core genes and lineage‐specific changes at peripheral genes.
Conclusions
5
This study shows that one of the most C_3_‐like members of subgenus Tillandsia, T. leiboldiana, is not fully C_3_, but a weak, facultative CAM plant. In contrast, T. vanhyningii exhibited strong physiological and transcriptomic drought tolerance consistent with its constitutive CAM phenotype and epiphytic lifestyle. Although these closely related species accumulated organic acids overnight under water limitation, they deployed largely distinct transcriptional architectures, particularly outside the core CAM (de)carboxylation module. This suggests that the evolutionary trajectory linking facultative and constitutive CAM involves broad transcriptional rewiring at peripheral parts of the networks, not simply the quantitative upregulation of a conserved gene set. The remarkable phenotypic and transcriptomic diversity within the subgenus Tillandsia highlights its potential as an outstanding system for exploring the CAM continuum and the evolutionary flexibility underlying this key innovation.
Author Contributions
This study was conceived by C.L., J.H., M.D.L.H., O.P. and C.G.‐C. Accessions were sampled by W.T. The greenhouse experiment and sampling were conducted by M.D.L.H. and M.H.J.B. Laboratory work including TA measurements and RNA extraction and library prep was conducted by M.D.L.H. and M.H.J.B. under the guidance of O.P. Analyses were performed by C.G.‐C., M.D.L.H. and S.S. Gas exchange measurements were done by G.B. with input from W.W. The manuscript was written by C.G.‐C. with initial contributions from S.S., was amended by O.P. and was finally read and approved by all coauthors.
Funding
This research was funded by the Austrian Science Fund (FWF) (grant DOI 10.55776/W1225 to a faculty team including C.L. and O.P., grant DOI 10.55776/P35275 to O.P.) and by the bilateral PRCI ANR‐FWF RadiaSpe project (ANR‐23‐CE02‐0032, FWF DOI: 10.55776/I6765). Additionally, this work was funded by the professorship start‐up grant of Christian Lexer at the University of Vienna BE772002. For open access purposes, the author has applied a CC BY public copyright licence to any author accepted manuscript version arising from this submission.
Conflicts of Interest
The authors declare no conflicts of interest.
Supporting information
Figure S1: Diagnostic plots for linear mixed‐effects models constructed in each species to assess the effect of time and watering regime on titratable acidity measurements. Figure S2: Mapping rates across 24 samples of different species and watering regimes. Figure S3: Principal component analysis of TMM‐normalised read counts in log(CPM). Figure S4: Overview of the overlap between DE analyses across time points and across watering regimes. Figure S5: Gene trees of PPC1 (A) and PPC2 (B). Figure S6: Gene trees of PPCK with the isozymes PPCK1 and PPCK2. Figure S7: Expression of MDH orthologs. Figure S8: Expression of NAD‐ME. Figure S9: Expression of PEPCK. Figure S10: Expression of PPDK and PPDK regulatory protein. Figure S11: Expression of carbonic anhydrases. Figure S12: Transcript abundance in ‘transcripts per million’ (TPM) for a subset of single‐copy orthogroups related to drought response in T. vanhyningii and T. leiboldiana.
Table S1: Sampling overview and summary of RNA extractions for 24 samples. Table S2: Full list of CAM‐related orthogroups obtained from Yardeni et al. (2021) based on A. comosus (pineapple) genes. Table S3: List of 61 CAM‐related orthogroups selected for expression visualisation in Figure 4. Table S4: List of 31 selected orthogroups related to water deprivation and drought response for expression visualisation in Figure S7. Table S5: Titratable acidity measurements for Tillandsia vanhyningii and Tillandsia leiboldiana at four time points spread over 6‐h intervals. Table S6: Statistical models and results for titratable acidity measurements. Table S7: Summary of per‐gene read counting with FeatureCounts. Table S8: Overlap of DE orthogroups between species, watering regimes and time points. Table S9: List of significantly enriched GO terms in each intersection of the Venn diagram showing overlap of DE genes between time points. Table S10: Description of CAM genes in sections of the Venn diagram A (day vs. night). Table S11: Description of CAM genes in sections of the Venn diagram B (water‐limited vs. standard watering conditions). Table S12: List of significantly enriched GO terms in each intersection of the Venn diagram showing overlap of DE genes between watering regimes. Table S13: Gas exchange measurements of T. leiboldiana under a progression from well‐watered to water‐limited conditions. Table S14: Gas exchange measurements of T. leiboldiana after a 14 days recovery following water‐withholding.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abraham, P. E. , H. Yin , A. M. Borland , et al. 2016. “Transcript, Protein and Metabolite Temporal Dynamics in the CAM Plant Agave.” Nature Plants 2, no. 12: 16178.27869799 10.1038/nplants.2016.178 · doi ↗ · pubmed ↗
- 2Alexa, A. , J. Rahnenführer , and T. Lengauer . 2006. “Improved Scoring of Functional Groups From Gene Expression Data by Decorrelating GO Graph Structure.” Bioinformatics 22, no. 13: 1600–1607.16606683 10.1093/bioinformatics/btl 140 · doi ↗ · pubmed ↗
- 3Arumingtyas, E. L. , E. S. Savitri , and R. D. Purwoningrahayu . 2013. “Protein Profiles and Dehydrin Accumulation in Some Soybean Varieties (Glycine max L. Merr) in Drought Stress Conditions.” American Journal of Plant Sciences 4, no. 1: 134–141.
- 4Assmann, S. M. , and T. Jegla . 2016. “Guard Cell Sensory Systems: Recent Insights on Stomatal Responses to Light, Abscisic Acid, and CO 2 .” Current Opinion in Plant Biology 33: 157–167.27518594 10.1016/j.pbi.2016.07.003 · doi ↗ · pubmed ↗
- 5Bachmann, A. , B. Hause , H. Maucher , et al. 2002. “Jasmonate‐Induced Lipid Peroxidation in Barley Leaves Initiated by Distinct 13‐LOX Forms of Chloroplasts.” Biological Chemistry 383, no. 10: 1645–1657.12452441 10.1515/BC.2002.185 · doi ↗ · pubmed ↗
- 6Barfuss, M. H. J. , W. Till , E. M. C. Leme , et al. 2016. “Taxonomic Revision of Bromeliaceae Subfam. Tillandsioideae Based on a Multi‐Locus DNA Sequence Phylogeny and Morphology.” Phytotaxa 279, no. 1: 1–97.
- 7Barghi, N. , J. Hermisson , and C. Schlötterer . 2020. “Polygenic Adaptation: A Unifying Framework to Understand Positive Selection.” Nature Reviews. Genetics 21, no. 12: 769–781.10.1038/s 41576-020-0250-z 32601318 · doi ↗ · pubmed ↗
- 8Barghi, N. , R. Tobler , V. Nolte , et al. 2019. “Genetic Redundancy Fuels Polygenic Adaptation in Drosophila .” P Lo S Biology 17, no. 2: e 3000128.30716062 10.1371/journal.pbio.3000128 PMC 6375663 · doi ↗ · pubmed ↗