Comparative transcriptome analysis of emerging young and mature leaves of Bienertia sinuspersici, a single-cell C4 plant

Plant material

Bienertia sinuspersici plants were maintained in 10-gallon (37.85-liter) citrus pots in growth chambers under a 14-hour light/10-hour dark cycle with a stepwise light regime increasing to 525 PPFM at full light and an 18 °C (dark) to 35 °C (light) temperature regime. Plants were watered once a week and were fertilized with Peters 20-21-5 in between waterings. Within two hours after light initiation, whole, fully expanded, mature leaves and newly emerging young leaves (approximately 0.2 mm—File S1) were harvested from three 9-month-old vegetative stage plants, combined as a pooled sample, and immediately flash-frozen in liquid nitrogen. The timing of sample harvesting coincided with the fully photosynthetic stage of the leaf. Samples were pooled from three plants due to the limited tissue available for RNAseq. Flash-frozen leaf tissue was ground into a fine powder with a liquid nitrogen-cooled mortar and pestle. Approximately 100 mg of frozen powder was transferred to a liquid nitrogen-frozen 2 mL Eppendorf tube and stored at −80 °C until RNA was extracted. Portions of the text used in this section were previously published as part of a preprint (Sharpe et al., 2023).

RNA extraction

Total RNA was extracted using an acid guanidinium thiocyanate phenol chloroform extraction method similar to that described previously (Chomczynski & Sacchi, 1987). Briefly, one mL of 0.8 M guanidinium thiocyanate, 0.4 M ammonium thiocyanate, 0.1 M sodium acetate pH 5.0, 5% w/v glycerol, and 38% v/v water saturated phenol were added to approximately 100 mg powdered tissue, shaken to evenly mix the sample, and incubated at room temperature for 5 min. 200 µL chloroform was added and shaken vigorously until the entire sample became uniformly cloudy before incubation at room temperature for 3 min. Samples were then centrifuged at 17,000× g at 4 °C for 15 min, and the aqueous phase was removed and transferred to a clean 1.5 mL Eppendorf tube. 600 µL 2-propanol was added, rocked 5 to 6 times, and incubated at room temperature for 10 min. Samples were centrifuged 17,000× g at 4 °C for 10 min, and the supernatant was discarded. one mL 75% DEPC-treated water mixed with ethanol was added, after which the samples were vortexed for 10 s and then centrifuged at 9,500× g at 4 °C for 5 min. The ethanol was removed by pipetting, and the pellets were allowed to dry completely. Pellets were suspended in RNase-free water and incubated at 37 °C with RNase-free DNase I (Thermo Scientific) for 30 min; the DNase I was inactivated at 65 °C for 10 min. 450 µL buffer RLC from the Qiagen (Valencia, CA, USA) RNeasy Plant Mini Kit was added to the digestion, processed in accordance with the manufacturer’s recommendations, and eluted in 50 µL RNase free water. Extracted RNA was quality checked either with the Bio-Rad (Hercules, CA, USA) Experion system using the Experion RNA High Sens Analysis kit or the Agilent (Santa Clara, CA, USA) 2100 Bioanalyzer system using the RNA Nano Chip and Plant RNA Nano Assay Class.

Illumina sequencing

cDNA and final sequencing library molecules were generated with Illumina’s TruSeq RNA Sample Preparation v2 kit and instructions with minor modifications. Modifications to the published protocol include a decrease in the mRNA fragmentation incubation time from 8 min to 30 s to create the final library proper molecule size range. Additionally, A7Biosciences’ (Woburn, MA, USA) DNA SizeSelector-I bead-based size selection system was utilized to target final library molecules for a mean size of 450 base pairs. All libraries were then quantified on a Life Technologies (Carlsbad, CA, USA) Qubit Fluorometer and qualified on an Agilent (Santa Clara, CA, USA) 2100 Bioanalyzer (Dr. Jeff Landgraf, personal communication). The Illumina HiSeq 2000 sequencing platform was used to sequence 2x100 PE reads from the cDNA libraries generated from the above RNA extractions at Michigan State University’s Research Technology Support Facility.

454 sequencing

cDNA libraries were constructed from the RNA extractions using the SMARTer™ PCR cDNA Synthesis Kit from ClonTech (Mountain View, CA, USA) according to the manufacturer’s instructions. cDNA quality and size distribution were verified via 1% TAE gels and the Bio-Rad (Hercules, CA, USA) Experion system. cDNA libraries were processed to attach the Rapid Library Multiplex Identification (RL MID) Adapters according to the manufacturer’s protocol. Libraries were then quality checked for size distribution with Agilent’s (Santa Clara, CA, USA) 2100 Bioanalyzer, quantified via fluorometry, pooled, and then sequenced on Roche Applied Science’s (Indianapolis, IN, USA) Genome Sequencer FLX System with GS FLX Titanium technology.

Data availability

All sequencing data was submitted to the NCBI SRA and was assigned accession number PRJNA340188.

Sequencing data QC and assembly

Sequence read information from Roche’s GS FLX Standard Flowgram Format (sff) files included 70,867 reads for the mature dataset and 54,462 reads for the emergent dataset; read information from Illumina HiSeq 2000 2x100 PE fastq files included 178,716,218 reads for the mature dataset and 218,726,388 reads for the emergent dataset. All developmental read datasets were processed with the CLC Create Sequencing QC Report tool to assess read quality. The CLC Trim Sequence process was used to trim the 454 read datasets for a Phred value of 15; the Illumina reads were trimmed for a Phred score of 30, and the first twelve 5′bases were removed due to GC ratio variability. All read datasets were trimmed of ambiguous bases. Illumina reads were then processed through the CLC Merge Overlapping Pairs tool. In the absence of an available reference genome for Bienertia sinuspersici, we performed a de novo assembly using CLC Genomics Workbench v6 and according to methods described previously (Hewitt, Ghogare & Dhingra, 2020; Sharpe et al., 2020; Hewitt et al., 2021). Trimmed reads used for assembly were mapped back to the assembled contigs, mapped reads were used to update the contigs, and contigs with no mapped reads were ignored. Consensus contig sequences were filtered to include only those with a consensus length of >200 and a depth of coverage of >5. The resulting dataset was exported as a multi-fasta file. The individual mature and emergent leaf read datasets, from the original non-trimmed reads, were mapped back to the assembled contigs to generate individual developmental sample read counts for each contig; read counts were then normalized with the reads per kilobase per million reads (RPKM) method (Mortazavi et al., 2008).

Functional annotation

The transcriptome fasta file produced from the assembly was imported into the OmicsBox Functional Annotation Module (Biobam Bioinformatics S.L., Valencia, Spain), following which contig sequences were identified by a blastx alignment against the NCBI ‘Viridiplantae’ database, gene ontology (GO) mapping was performed, and annotations were assigned using the Blast2GO feature with default parameters (Götz et al., 2008).

Differential expression analysis

Pairwise differential expression analyses were conducted for the emergent versus mature samples via the OmicsBox Transcriptomics Module using NOISeq-sim to compare the two developmental stages. NOISeq-sim infers significant differential expression without experimental replicates. For absent replicates, NOISeq-sim uses a multinomial distribution to model technical replicate read counts (Tarazona et al., 2011; Tarazona et al., 2013; Tarazona et al., 2015). Our research group has employed this method successfully in previous studies where it was not possible to obtain more than a single replicate per sample (Hewitt et al., 2021; Hewitt et al., 2023). Bienertia sinuspersici is difficult to sample in the wild, as it grows primarily in areas around the Persian Gulf. Since it is adapted to hot, dry, and high salt environments, it isn’t easy to grow, even in climate-controlled chambers. Default parameters were used to simulate five replications with a set variability of 0.02 in each replication. Genes with a NOIseq probability greater than 0.9 and a —log2 fold change expression— value greater than 1.0 for at least one sample group were considered differentially expressed. A complete list of DE genes (along with the non-DE genes) and their corresponding functional annotations and expression values can be found in File S2.

GO enrichment analysis

GO enrichment analysis using Fisher’s exact test was conducted in the OmicsBox Functional Annotation Module to identify the cellular components, molecular functions, and biological processes that were enriched in each of the two developmental stages. Lists of the differentially expressed, functionally annotated genes were generated for the emergent and mature Bienertia. These lists served as the treatment datasets for enrichment analyses, and the master annotated transcriptome was used as the reference dataset. In addition to conducting a full GO enrichment analysis, a second analysis was run using the GO-Slim feature to reduce the number of GO terms present in the annotated reference transcriptome to include only general functions and processes displaying the greatest enrichment. For both analyses, an FDR-corrected p-value of 0.01 was used as the cutoff for determining statistical significance and to reduce the GO assignations to the most specific terms. Complete GO enrichment analysis and GO-slim analysis results can be found in File S3.

Gene expression and gene ontology enrichment analyses

De novo assembly yielded total of 116,257 contigs assembled from 141,504,502 trimmed reads with an N50 of 792 bases. Filtering for contigs with consensus length >200 and >5x coverage yielded 72,524 transcripts. Top blast hit descriptions were assigned to 26,007 contigs, and a subset of 24,603 contigs were fully annotated with corresponding gene ontologies and interpro IDs (File S2).

Expression analysis identified 3,364 DE genes, 1,761 in the emergent tissue and 1,603 genes in the mature tissue with high probability of differential expression. Of these, 1,062 and 661 DE genes were matched with functional annotations in the emergent and mature tissues, respectively (File S2). Emergent leaves exhibited slightly higher median and mean global and DE gene expression levels than mature leaves, as reflected in the log10-transformed RPKM values (Fig. 1). Additionally, a broader range of expression values was observed in the emergent leaves compared to the mature leaves, suggesting more diverse transcriptional activity in the former. This observation may reflect a broader range of cellular activities required for early leaf development and cellular differentiation.

The assignment of GO terms to annotated datasets enabled the identification and classification of biological processes, molecular functions, and cellular components that were overrepresented in each of the treatment datasets versus the entire transcriptome. A total of 110 GOs associated with biological processes (bpGOs), 43 GOs associated with molecular functions (mfGOs), and 21 GOs associated with cellular components (ccGOs) were uniquely overrepresented in the emergent tissues; 48 bpGOs, 23 mfGOs, and six ccGOs were uniquely overrepresented in the mature tissues; and 27 bpGOs, 11 mfGOs, and 10 ccGOs were overrepresented in both emergent and mature tissues (File S3). To further simplify the dataset and identify the more inclusive ontologies representative of each developmental stage, we used the OmicsBox GO-slim tool to create a reduced list of enriched ontologies. The simplified ontology list was comprised of 10 bpGOs, four mfGOs, and four ccGOs uniquely enriched in the emergent tissue; two bpGOs, two mfGO, and one ccGO uniquely enriched in the mature tissue; and 20 bpGOs, seven mfGOs, and eight ccGOs that were enriched in both emergent and mature tissues (File S3).

GO-slim results provided a broad overview of enriched functions, guiding the analysis by highlighting key pathways. More specific GO terms (from the comprehensive enrichment analysis) corresponding to the GO-slims offer detailed insight into the enriched processes at each developmental stage (Table 1). In emergent tissues, unique enrichment was observed in terms related to cell cycle, growth, and differentiation, including specific GOs such as ‘cyclin-dependent protein serine/threonine kinase activity’, ‘regulation of G2/M transition of mitotic cell cycle’, ‘regulation of cell division’, ‘regulation of cell growth/size’, ‘asymmetric cell division’, and ‘regulation of cell population proliferation’. Terms associated with cytoskeleton and cytoskeletal motor activities were also uniquely enriched in the emergent stage, with specific GOs including ‘cortical microtubule’, ‘microtubule-based movement’, ‘spindle assembly’, and ‘phragmoplast’. These results highlight the stage-specific importance of processes essential for cell division, structural establishment, and spatial organization, which are necessary for preparing emergent tissues for later specialization.

In contrast, mature tissues exhibited unique enrichment of GOs associated with photosynthetic processes, including specific terms such as ‘chloroplast thylakoid membrane’, ‘chlorophyll binding’, ‘photosystem II, light harvesting’, and ‘photosynthetic electron transport chain’. Terms related to homeostasis and oxygen-binding processes, such as ‘oxygen binding’ and ‘peroxidase activity’, were also uniquely enriched, reflecting the tissue’s transition to specialized photosynthetic functions and adaptations for maintaining internal stability. Additionally, carbohydrate-related terms, such as ‘carbohydrate binding’ and ‘sucrose-phosphate synthase activity’, highlight the importance of resource management and structural integrity in photosynthetically active tissue.

GO terms that were highly enriched at both developmental stages included terms associated with transcription-mediated response to stimuli—e.g., ‘DNA binding’, ‘DNA-binding transcription factor activity’, and ‘regulation of DNA-templated transcription’—and cellular component organization—e.g., ‘system development’, ‘cytoplasmic vesicle’, ‘vacuole’, ‘mitochondria’, and ‘chloroplast stroma’. The shared enrichment of these GOs highlights core processes critical for transcriptional regulation and organelle organization at both stages of development.

Analysis of key DE genes corresponding to enriched pathways provides further insights into the developmental processes driving SCC₄ functionality in Bienertia. The following sections examine key genes associated with the enriched GO terms, highlighting the roles these genes play as tissues transition from growth and differentiation in the emergent stage to functional specialization in the mature stage.

Genes associated with the cell cycle and growth

In early leaf development, meristematic cells actively divide; however, as leaves mature, cell proliferation becomes confined to the leaf base while cells elongate and expand (Gonzalez, Vanhaeren & Inzé, 2012; Kalve, De Vos & Beemster, 2014). Underpinning these developmental processes, the cell cycle is regulated by cyclins (CYCs), cyclin-dependent kinases (CDKs), aurora kinases, condensin complexes, and mitotic spindle checkpoints, which collectively govern the transitions between cell cycle phases (Collette et al., 2011; Willems et al., 2018; Shimotohno et al., 2021; Li et al., 2024a; Li et al., 2024b).

Aligning with this established information, in Bienertia, the largest group of enriched GO terms identified in the emergent tissue were associated with cell cycle regulation, cell growth, and cellular differentiation. Analysis of the DE genes corresponding to these enriched GOs revealed differential expression of genes corresponding to various components of the cell cycle, particularly genes associated with cellular duplication. Transcripts corresponding to CDKB1-1, along with CYCA1, CYCA2, CYCA3-1, CYCD2, CYCD3, CYCA1-like, and CYCS13-6 were among those that were significantly upregulated in the emergent tissue (Table 2); all of these are known to accumulate during the mitotic synthesis phase (S), the second gap (G2) phase, and the mitosis (M) phase (Boudolf et al., 2004; Boudolf et al., 2009; Inzé & De Veylder, 2006). Only one cyclin-associated transcript, CDKF4, was differentially expressed in mature tissue. Along with CYC/CDK, the emergent tissue displayed significantly increased expression of transcripts corresponding to serine-threonine protein kinases AURORA1, AURORA2, AURORA3 (AUR1, AUR2, AUR3); CONDENSIN COMPLEX SUBUNIT 1, 2, and 3 (SMC1, SMC2, and SMC3); and mitotic spindle checkpoint-associated MITOTIC ARREST DEFICIENCY 2 and 3 (MAD2 and MAD3). AUR and SMC are important mediators of cellular mitosis, with crucial roles in G2/M transition, chromosome binding, and kinetochore separation (Collette et al., 2011; Willems et al., 2018).

In addition to core cell cycle regulators, the developmental patterning of leaf tissue can be monitored via the expression of leucine-rich repeat receptor-like protein kinases (LRR-RLK) family members. These proteins detect and transduce signals to initiate responses critical for the development of shoot organs. One clade of the LRR-RLK family, including ERECTA and ERECTA-like (ERL), localizes and functions in the shoot apical meristem (SAM) and organ primordia. Previous work has demonstrated that ERECTA mRNA is expressed at low levels in the SAM, with expression increasing in developing vegetative organs and decreasing in mature organs (Yokoyama et al., 1998). Consistent with this, Bienertia emergent leaf tissue displayed significantly higher expression levels of transcripts corresponding to ERECTA and ERL compared to mature tissues (Table 2). Not all LRR-RLK genes are developmentally expressed to promote actively dividing cells in emergent tissues; some, such as PHLOEM INTERCALATED WITH XYLEM 2-like (PXL2-like) and wall-associated receptor kinase CLAVATA1 (CLV1), act in maturing tissues to suppress cell division or generate function-specific cell types in expanding tissues. PXL2-like has been implicated in the development of phloem and xylem (Etchells et al., 2013), while CLV1 facilitates phyllotaxis formation in leaf primordia and suppresses undifferentiated cells at the shoot meristem, committing them to organ development (Clark, Running & Meyerowitz, 1993). These genes, significantly upregulated in mature tissues, suggest an active role in vascular differentiation and morphological organization during the later stages of leaf development. Beyond the LRR-RLK family, wall-associated kinases (WAKLs) contribute to the development of leaf structure and morphology during development. Several WAKL genes (WAKL2, WAKL4, WAKL8, and WAKL9), which are known to mediate cellular signaling and wall dynamics, were highly expressed in mature tissues (Verica et al., 2003; Sharma et al., 2020). Together, these findings support the structural and physiological adaptations necessary for functional specialization in mature leaves.

Genes involved in cytoskeletal and structural organization

Cytoskeletal dynamics are central to the coordination of cell division and structural development, providing stability and facilitating key processes such as mitotic spindle formation, organelle positioning, and cytoplasmic streaming (Smertenko et al., 2018). The cytoskeleton consists of actin filaments and microtubules, each serving distinct functions: actin primarily drives cytoplasmic streaming, while microtubules support cell wall assembly and intracellular transport. Microtubule-associated proteins (MAPs), including motor MAPs like kinesins and non-motor MAPs that stabilize and regulate microtubule dynamics, are crucial for these functions (Parrotta, Cresti & Cai, 2014; Lee, Qiu & Liu, 2015). In Bienertia, cytoskeletal components are particularly relevant to the development of dimorphic chloroplasts, as the chloroplasts of Bienertia species have been shown to interact with microtubules and actin as the chlorenchyma cell develops (Voznesenskaya et al., 2005; Chuong, Franceschi & Edwards, 2006) (Table 1).

Nearly all differentially expressed cytoskeleton-associated genes were more highly expressed in emergent leaf tissues, highlighting the critical role of cytoskeletal components in cell division, structural establishment, and intracellular organization during early leaf development. β-ACTIN (ACTB), a key component of actin filaments, was highly expressed at both developmental stages with significant differential expression in emergent tissues. The comparatively high expression of ACTB in early development may reflect its importance for organelle positioning and the establishment of foundations for further compartmentalization. Notable among the motor MAPs expressed in emergent tissues are kinesin-like proteins (KIFs), such as KIF18B, KIF4A, and KIF22, which are involved in spindle dynamics, chromosome alignment, and intracellular transport, processes integral to cell division and structural organization (Lee, Qiu & Liu, 2015). Additionally, non-motor MAP genes such as MAP65-3 and MAP65-5 contribute to microtubule bundling and stabilization during mitosis, facilitating cytokinesis and cell wall assembly (Parrotta, Cresti & Cai, 2014). α-TUBULIN (TUBA) and β-TUBULIN3 (TUBB3) were similarly expressed, highlighting their involvement in spindle formation and intracellular movement during cell proliferation. Together, these findings highlight the important role of microtubule dynamics, not only in establishing cellular architecture but also in enabling the structural and functional transitions unique to Bienertia.

Genes associated with photosynthetic development

Photosynthetic development in Bienertia is governed by transcriptional programs that drive chloroplast biogenesis, maturation, and the establishment of dimorphic chloroplast functionality required for SCC₄ carbon fixation. Most photosynthesis and carbon fixation-related genes identified in this study showed differential expression in mature tissues (Table 1), consistent with the enriched GO category for photosynthesis, while genes in emergent tissues appeared to establish foundational processes.

In mature tissues, genes encoding light-harvesting complex proteins, such as LIGHT HARVESTING CHLOROPHYLL-BINDING (LCHB), LCHB-like, and LCHB6, facilitate energy capture and transfer to the photosystems, critical for sustaining photosynthesis. Genes encoding subunits of PHOTOSYSTEM I (PSI)—including PSA-C, PSA-E, PSA-H, PSA-I, and PSA-K—were more highly expressed in mature tissues, highlighting their roles in electron transfer and energy conversion at this stage (Jensen et al., 2000; Jensen et al., 2007). Regarding PHOTOSYSTEM II (PSII), the D1 PROTEIN encoding gene (PSII D1), a core component involved in water splitting and oxygen evolution, and PSB-Y, a stabilizer of the PSII complex, were also strongly expressed. The mature tissue-specific expression of IRON SUPEROXIDE DISMUTASE (FeSOD) likely contributes to redox balance within chloroplasts, mitigating oxidative stress generated during high rates of photosynthesis. Regulatory and biosynthetic processes supporting chloroplast functionality were evident in mature tissues. RNA POLYMERASE SIGMA FACTOR D (RPO-D) and SIGMA FACTOR (SIG) were expressed, potentially regulating chloroplast gene expression. MAGNESIUM-PROTOPORPHYRIN IX MONOMETHYL ESTER CYCLOHYDROLASE (MgPME cyclase), which supports chlorophyll biosynthesis, was also highly expressed in mature tissues, potentially reflecting an increased demand for chlorophyll to sustain higher rates of photosynthesis in mature leaves (Takahashi et al., 2014; Kong et al., 2016).

While photosynthetic gene expression in emergent tissues was generally lower than in mature tissues, several genes essential for early chloroplast biogenesis and organization that were expressed at the emergent stage included LIGHT-INDUCED PROTEIN 3 (LIL3), which supports chlorophyll production by stabilizing biosynthetic intermediates (Takahashi et al., 2014). Structural assembly and functional development of emergent chloroplasts appears to be further supported by OUTER ENVELOPE PROTEIN 80 (OEP80) and TRANSLOCON OF THE OUTER CHLOROPLAST MEMBRANE 75 (TOC75), which mediate the import of nuclear-encoded proteins into chloroplasts (Yoshimura et al., 2024). These genes may play critical roles in establishing the distinct structural and functional identities of central and peripheral chloroplasts in Bienertia.

With regards to genes associated with carbon fixation, high expression of CA was observed in mature tissues, reflecting its role in supporting efficient carbon fixation within peripheral chloroplasts as SCC₄ photosynthesis becomes established. Additionally, PEPC was expressed at both stages, with one transcript more highly expressed in emergent tissues and another in mature tissues. This suggests that the differential expression of PEPC transcripts may play a role in supporting the establishment and specialization of distinct chloroplast functions during the transition from C₃ to SCC₄ photosynthesis. Specifically, these expression patterns could reflect the progressive development of the spatially separated roles of peripheral and central chloroplasts in the SCC₄ pathway. Finally, the elevated expression of two transcripts corresponding to the RUBISCO small subunit in mature tissues reflects its increased role in central chloroplasts during the mature stage, when carbon fixation becomes fully operational. This pattern suggests a developmental shift in chloroplast function to support the heightened demands of SCC₄ photosynthesis.

Differentially expressed transcription factors

Among the GOs enriched at both emergent and mature stages were several associated with transcriptional regulation. Given the central role of transcription factors (TFs) in orchestrating developmental processes, stage-specific expression of TFs in Bienertia was explored (Table 3).

In emergent tissues, a higher diversity and abundance of transcription factor-associated genes were observed, including genes encoding auxin response factor (ARF), basic helix-loop-helix (bHLH), homeodomain (HD), MYB, squamosa promoter-binding protein-like (SPL), and zinc finger (ZF) family proteins. These TFs likely regulate gene expression programs critical for cell proliferation, differentiation, and early leaf anatomical development. For example, ARF and MYB transcription factor genes, which are implicated in hormonal signaling pathways such as auxin and gibberellin signaling (Shin et al., 2007; Zhao et al., 2021), were highly expressed in emergent tissues, suggesting a role in the coordination of hormone responses necessary for growth and the establishment of leaf structures. Similarly, NAC and SPL transcription factor genes, known to mediate stress responses and metabolic regulation (Filichkin et al., 2018; Hernandez, Goswami & Sanan-Mishra, 2020; Min et al., 2022), likely contribute to early adaptation mechanisms ensuring survival and proper development during this critical stage, while BEL1-like, a homeodomain-containing TF, functions in the shoot apex to regulate meristem identity and facilitate the correct development of shoot architecture (Kumar et al., 2007).

The heightened expression of bHLH transcription factors, such as TRANSCRIPTION FACTOR SPEECHLESS-LIKE, at the emergent stage, aligns with GO terms related to cell cycle regulation and cytoskeletal dynamics, pointing to their involvement in stomatal developmental progression and structural organization (Chen, Wu & Hou, 2020). Additionally, GATA TF genes, with known roles in regulating chloroplast biogenesis (Chiang et al., 2012; Hudson et al., 2013), may lay the groundwork during the emergent stage for the transition to photosynthetic functionality.

In mature tissues, fewer TF families were differentially expressed, but a distinct enrichment of highly expressed ethylene-responsive transcription factor genes (ERFs) was observed. These ERFs, including ERF1A, ERF2, and DEHYDRATION-RESPONSIVE ELEMENT BINDING PROTEIN, are commonly associated with stress responses, senescence, and hormone-mediated cross-talk (Nakano et al., 2006). The prominence of ERF genes suggests their regulatory roles in maintaining cellular homeostasis and responding to environmental stimuli during the mature photosynthetic stage. Additionally, the mature stage showed expression of CONSTANS (CO) and CONSTANS-like (COL), zinc finger TFs linked to photoperiodic regulation (Romero et al., 2024), which aligns with the enriched GO terms for photosynthesis and light response.

Framework for SCC₄ and chloroplast differentiation

The findings of our transcriptomic analysis provide a framework for understanding how stage-specific gene expression contributes to the structural and functional transitions underlying the shift from C₃ photosynthesis in young Bienertia leaves to the specialized SCC₄ photosynthetic mechanism in mature leaves (Fig. 2).

In emergent tissues, high expression of cell cycle regulators such as cyclins, cyclin-dependent kinases, and aurora kinases, along with cytoskeletal components like tubulins, kinesins, and microtubule-associated proteins, indicates the occurrence of active cell division and cytoplasmic reorganization. These processes are likely critical for establishing the cellular architecture required for chloroplast partitioning. Genes such as TOC and OEP80, involved in protein import into chloroplasts, suggest that foundational steps in chloroplast biogenesis, and perhaps differentiation, are initiated during this stage. Additionally, expression of BEL1-like regulates meristem identity, suggesting broader developmental signals coordinating shoot architecture and cellular differentiation. These findings suggest that early developmental processes active in emergent tissues lay the groundwork for the physical and functional separation of chloroplast types.

As leaves mature, the transcriptional profile shifts toward processes supporting SCC₄ functionality. Photosystem-related genes, including components of PSI, PSII, and light-harvesting complexes, are highly expressed, reflecting their roles in optimizing light capture and energy conversion. This is complemented by the expression of regulatory genes like SIG and RPO-D, essential for chloroplast-specific gene expression, as well as MgPME cyclase, which supports chlorophyll production. The upregulation of FeSOD indicates that oxidative stress management becomes particularly important in mature tissues, where high rates of photosynthesis generate reactive oxygen species.

Carbon fixation-associated enzymes also exhibit stage-specific expression patterns that shed light on the C₃–SCC₄ transition. The differential expression of PEPC transcripts in both emergent and mature tissues suggests functional specialization of this gene during development, with distinct isoforms potentially aligning with early metabolic priming and later bicarbonate fixation in peripheral chloroplasts. Similarly, elevated expression of CA and the RUBISCO small subunit in mature tissues reflects their roles in establishing the spatial and functional compartmentalization of SCC₄ photosynthesis, with peripheral and central chloroplasts specializing in initial and secondary CO₂ fixation, respectively.

Overall, this model highlights how the emergent stage establishes the cellular framework and initiates key molecular processes, while the mature stage refines and specializes these systems for efficient light harvesting and carbon fixation.