Tender leaves of R. sericea from the same individual autumn shoots (Figures 1A, B), distinct from the spring leaves, were collected in triplicate on September 28, 2023, at 2,518 meters, on Erlang Mountain (Luding County, Sichuan Province, China). The leaves were immediately preserved by rapid freezing in liquid nitrogen to maintain cellular integrity.

The nucleus isolation protocol was adapted from a previously published study (Conde et al., 2021). Frozen tissue was carefully ground into small pieces in liquid nitrogen using a mortar and pestle, and then transferred to a gentleMACS M tube (MiltenyiBiotec) containing 5 mL of Nuclei Isolation Buffer (NIB) (5% Dextran T40, 0.4M sucrose, 10 mM MgCl2, 1 mM DTT, 0.1% Triton X-100, 2U/μl Protector RNase Inhibitor, 100 mM Tris-HCl pH 7.4). The resulting suspension was filtered through a 70-µm mesh strainer and centrifuged at 300g for 1 minute. The pellet was resuspended in 500 µL NIB, filtered again through a 40-µm strainer, and subsequently resuspended in Wash Buffer (10 mM PBS, 1% BSA, 2U/μl Protector RNase Inhibitor). To minimize contamination from organelles such as chloroplasts and to remove cellular debris that could potentially clog microfluidic chips, the nuclei were sorted by gating on the DAPI peaks using a BD FACS Aria III (BD Biosciences). The sorted nuclei were stained with 10 µM DAPI for 5 minutes to assess their integrity and concentration.

The snRNA-seq libraries were constructed using the Chromium Next GEM Single Cell 3’ GEM, Library & Gel Bead Kit v.3.1, following the instructions provided in the user manual (Figure 1C). Additionally, three separate preparations of single-nucleus suspensions, derived from tender leaves of the R. sericea shrub, were sequenced to represent three biological replicates.

After completing the library preparation, quality control was performed. Initial quantification was done using the Qubit 2.0 fluorometer, followed by assessment of the insert DNA size with the Agilent 2100 Bioanalyzer. Once the insert size met the expected criteria, the effective library concentration (2nM) was accurately quantified using qPCR to ensure the library’s quality. After passing quality control, the libraries were sequenced on the Illumina NovaSeq 6000 platform using 150 bp paired-end reads. The sequencing depth of each sample was at least 100G.

Before bioinformatics analysis, the raw sequencing data were assessed for quality using FastQC. The data were then subjected to several quality control and filtering steps to generate clean data for further analysis.

The Cell Ranger software v3.1.0 from 10× Genomics was utilized to align sequencing reads to the reference genome, assess cell numbers, and filter cells (Zheng et al., 2017). Initially, we performed read alignment using STAR (Dobin et al., 2013) to map the sequencing reads to the R. sericea reference genome (CNA0146539). The uniquely aligned reads were then quantified to determine UMI types and generate the UMI-barcode expression matrix. Subsequently, cellular and non-cellular barcodes were distinguished by analyzing variations in unique molecular identifiers (UMIs) counts, which are indicative of gene expression levels across individual barcodes, which were assigned as candidate cell identifiers.

To ensure the retention of high-quality cells, the following criteria were applied: cells with fewer than 20% mitochondrial genes, cells expressing more than 200 genes, total expression levels greater than 400, and predicted doublet scores below 0.25, using the Python package Scrublet (Wolock et al., 2019). This filtering process successfully eliminated low-quality cells as well as doublets and multiplets. After applying these quality control criteria, 31,796 high-quality single cells were retained for subsequent analyses.

The Seurat v3 software (Stuart et al., 2019) was used to perform cell clustering, dimensionality reduction and differential gene expression analysis. We applied “FindIntegrationAnchors” function and “IntegrateData” function to remove the batch effects between samples. The matrix was standardized using the LogNormalize method from the “Normalization” function. To identify HVGs, the “Variance Stabilizing Transformation” method was applied in the “Find Variable Features” function, and the top 2,000 HVGs were selected. Cell clustering was performed using the “Find Clusters” function, with a weighted graph-based method called Shared Nearest Neighbour (SNN). For cell visualization, t-distributed stochastic neighbor embedding (t-SNE) was applied in Seurat (Kobak and Berens, 2019), and the results were concurrently validated using uniform manifold approximation and projection (UMAP) software (Becht et al., 2018).

Orthologous gene alignments of published leaf-related marker genes from the model plant Arabidopsis thaliana were utilized to identify cell-type-specific markers in R. sericea. Specifically, the detailed sequences of these marker genes were retrieved from the Plant Cell Marker Database (http://www.tobaccodb.org/pcmdb/) (Jin et al., 2022). Additionally, data from PlantscRNAdb (http://ibi.zju.edu.cn/plantscrnadb/#/) (Chen et al., 2021) and scPlantDB (https://biobigdata.nju.edu.cn/scplantdb/home) (He et al., 2024) were integrated for further supplementation and cross-verification. Using the A. thaliana (TAIR10) marker genes as query sequences, homologous genes in R. sericea were identified using the BLAST v2.12.0 (Camacho et al., 2009). The top-scoring hits were selected and annotated as corresponding R. sericea cell-type-specific genes.

To assess the robustness of cell type annotation without additional experimental validation, we conducted a cross-species transcriptomic comparison by mapping R. sericea leaf snRNA-seq data to a published A. thaliana leaf snRNA-seq reference atlas. We used the Rosette S3 dataset from a recent high-resolution study (Guo et al., 2025), which includes well-characterized leaf cell populations comparable to our annotations. Putative orthologs were identified by BLAST, retaining best-hit gene pairs with the highest bit score. Highly variable genes were defined independently in both datasets using Seurat (Stuart et al., 2019). Cross-species similarity was evaluated using Spearman’s rank correlation coefficient and visualized as a heatmap in R.

The hdWGCNA package (Morabito et al., 2023) in R (https://smorabit.github.io/hdWGCNA/) was employed for gene co-expression analysis. A weighted adjacency matrix was constructed using an unsupervised hierarchical clustering method (Morabito et al., 2021, 2023), with the optimal soft threshold power (β) was set to be 6 to assess the scale-free topology. Modules were identified based on the following parameters: gene_select = “fraction,” fraction = 0.05, nearest-neighbors parameter (k) = 50, minModuleSize = 50, and max_shared = 10. Significant co-expression modules were defined when the Pearson correlation coefficient was greater than 0.40 and the p-value was below 0.05. Finally, the co-expression network of genes was visualized using Cytoscape software v3.9.1 (Shannon et al., 2003).

The clusterProfiler R package was employed to perform the enrichment analyses (Yu et al., 2012). GO enrichment analysis (Gene Ontology, 2008) was conducted on differentially expressed genes (DEGs) to identify their potential functional enrichments. Based on Module Eigengene Connectivity (kME), hub genes identified from the hdWGCNA analysis were selected for the GSEA (Subramanian et al., 2005). Reference gene sets, such as the target gene sets (G6) from the PlantGSAD database were utilized (Ma et al., 2022). The parameters for the enrichment analysis were configured as follows: pAdjustMethod = ‘BH’, pvalueCutoff = 1, and qvalueCutoff = 1.

Differentiation trajectories of all leaf cells and leaf epidermal cells were analyzed using the Monocle 2 v2.26.0 R package to investigate their pseudotime relationships (Qiu et al., 2017). The “dpFeature” function was used to identify genes that define specific biological processes. Clustering analysis was performed using the “reduce Dimension ()” function, with max_components = 2 and reduction_method = “DDRTree.” The differentiation trajectories were then inferred using the “orderCells” function with default parameters (Qiu et al., 2017). Visualizations of gene expression were generated using the specialized pseudotime trajectory function to track changes in gene expression throughout differentiation stages. This method was also applied to several marker genes relevant to developmental pseudotime. Branch-dependent differentially expressed genes were identified using the BEAM function.

Total RNA was isolated from fresh R. sericea leaves, and the full-length coding sequence (CDS) of the RO6G37307 gene was amplified through reverse transcription. The overexpression vector was constructed via homologous recombination. A constitutive 35S promoter-driven overexpression vector was utilized for cloning. After homologous recombination, 5 µL of the recombinant product was introduced into Escherichia coli, which was then plated on LB agar containing 100 mg/L spectinomycin for colony selection. Four single colonies were picked for PCR verification, and plasmids were extracted from those with correct sequencing results. Subsequently, the Agrobacterium-mediated floral-dip method was applied to introduce the overexpression construct into A. thaliana (Columbia ecotype, Col-0) at the bolting stage. Plants were dipped three times at 5–7 day intervals to improve transformation efficiency. Seeds collected from the T0 generation were sown on kanamycin-containing medium to select positive transformants.

After obtaining the resistant lines through selection for kanamycin resistance, total DNA and RNA were extracted from the candidate plants using the BGMG Fast Total RNA/DNA Co-Extraction Kit. The resistant lines were then subjected to PCR and RT-qPCR validation. Stable RO6G37307 overexpression lines were established in the T2 generation. The leaf surface trichome of the RO6G37307-OE lines were imaged in three biological replicates using a digital camera (Nikon D750) and scanning electron microscopy (SEM) at three positions (base, mid, apex) on the 3rd and 5th rosette leaves, and the 1st cauline leaf. Each leaf section measured 5 mm × 5 mm, and the trichome density was compared to that of the Col-0 ecotype.

Homologous amino acid sequences corresponding to the RO6G37307 protein were retrieved from the NCBI database and filtered according to the following thresholds: E-value ≤ 1 × 10^-5, sequence identity ≥ 50%, and query coverage ≥ 70%. A phylogenetic relationship was inferred using the neighbor-joining (NJ) method implemented in MEGA v12.0 under the p-distance model, with 5,000 bootstrap replicates to assess node support; all other parameters were retained at their default values (Kumar et al., 2024b). Signal peptide prediction was conducted using SignalP v6.0 under default settings (Teufel et al., 2022), while transmembrane domain topology was assessed via DeepTMHMM v1.0 with default configurations (Hallgren et al., 2022).

For experimental subcellular localization, the full-length CDS of the RO6G37307 gene was fused in-frame to the enhanced green fluorescent protein (EGFP) within the pBI121 vector. The recombinant construct was introduced into one-month-old Nicotiana benthamiana leaves through Agrobacterium-mediated transient transformation (strain GV3101 carrying the pSoup-p19 helper plasmid). The infiltrated leaves were cultured at room temperature in darkness for 36 to 48 h. Leaf segments were then excised and mounted on glass slides. To increase the visibility of the plasma membrane and cell wall interface, samples were treated with 1 M mannitol solution for approximately 10 minutes to induce plasmolysis. The intracellular distribution of the protein was observed by observing EGFP fluorescence using a confocal laser scanning microscope.

To establish a comprehensive single-cell atlas capturing the developmental dynamics of autumn leaves in R. sericea for reflecting climate-induced secondary growth. we collected tender leaves from autumn shoots (Figures 1A, B), isolated nuclei, and performed snRNA-seq using the 10× Genomics Chromium platform. Following the standard snRNA-seq workflow (Figure 1C), we generated transcriptomic profiles from 9,368, 10,434, and 12,156 single cells across three independent biological replicates, which were subjected to rigorous quality control. On average, approximately 1,100 genes and 1,900 UMIs were detected per cell (Supplementary Table S1). After filtering out low-quality cells and potential doublets, 9,365, 10,357, and 12,074 high-quality single-cell transcriptomes were retained, yielding a final dataset of 31,796 cells for downstream analysis (Supplementary Table S2).

For cell population identification, we performed batch effect correction and integrated the three replicates. This analysis resolved 22 distinct cell clusters, which were visualized by UMAP (Figure 1D) and validated with t-SNE (Supplementary Figure S1). The clustering results were highly reproducible among replicates, demonstrating minimal batch effects and robust data quality. Furthermore, differential gene expression analysis and inter-cluster Pearson’s correlation further confirmed the reliability and biological relevance of the clustering results (Figures 1E, F).

To obtain reliable annotations, we first performed GO enrichment analysis on each cell cluster, aiming to infer their potential cellular identities based on biological functions. For instance, Clusters 0, 1, 4, 6, and 12 were predominantly enriched for chloroplast-related and photosynthesis functions (Supplementary Figure S2), suggesting their potential identity as mesophyll cell. Subsequently, we examined the top 50 cluster-specific genes and compared them with previously reported marker genes in A. thaliana (Table 1; Supplementary Table S3, Supplementary Figure S3). Furthermore, we analyzed the spatial expression patterns of 52 specific marker genesacross 22 cell clusters to more intuitively assess gene-specific expression and aid in leaf celltype annotation. (Supplementary Figure S4) This approach enabled the preliminary annotation of nine major cell types: mesophyll, epidermis, meristematic, proliferating, phloem, phloem parenchyma, xylem, bundle sheath, and companion cells (Figures 2A, C, D; Supplementary Figure S5).

As the outermost interface between plants and the aerial environment, the leaf epidermis provides a more sensitive indicator of environmental changes. However, the epidermal cell populations remained insufficiently resolved based on the first round of annotation results, which limits our understanding of the plateau adaptation of R. sericea. Therefore, we reanalyzed this subset to achieve higher resolution. A total of 11,888 epidermal cells were reclustered into 23 subclusters, which were subsequently putatively annotated as epidermal pavement, guard, and trichome cells (Figures 2B, E, F).

Integrating the results from both rounds of annotation, we identified 11 distinct leaf cell types: mesophyll, meristematic, proliferating, phloem, phloem parenchyma, xylem, bundle sheath, companion, epidermal pavement, guard, and trichome cells. This refined classification highlights the cellular diversity of R. sericea leaves and serves as a framework for downstream functional analyses.

To evaluate the consistency of our cell type annotations, we compared the R. sericea leaf snRNA-seq dataset with a published A. thaliana leaf snRNA-seq reference atlas using a cross-species transcriptomic similarity analysis. Based on orthologous highly variable genes, average expression profiles were calculated for each annotated cell type, and absolute Spearman rank correlation coefficients (r_s) were computed between different cell types of two species. Overall, several biologically corresponding cell type pairs showed modest but statistically significant expression similarity across species (r_s ≥ 0.24, p < 0.01), including companion, bundle sheath, xylem, mesophyll, and epidermis-related cells. Although the strength of correlation varied among cell types, the overall correspondence pattern was broadly consistent with expectations derived from marker-based annotations. Notably, certain cell populations displayed partial overlap with multiple A. thaliana reference types or relatively weaker similarity signals (Supplementary Figure S6). The relatively modest absolute correlation values observed here highlight both species-specific transcriptional programs, evolutionary divergence, and developmental context between the two datasets. Together, these results provide additional, transcriptome-level support for the plausibility of the assigned cell identities, while acknowledging inherent limitations associated with cross-species comparisons.

Beyond marker gene annotation, functional enrichment provides an additional layer of evidence for cell identity. To examine the classification and explore biological specialization, we performed GO enrichment analysis on the DEGs from each cell type. The enriched GO terms aligned closely with expected functional roles, with results consistent with the annotation and informative for adaptive processes (Supplementary Figure S7).

Specifically, mesophyll cells were enriched for photosynthesis-related categories such as “photosynthetic process” (GO:0015979) and “chloroplast thylakoid membrane” (GO:0009535). Meristematic and proliferating cells were associated with GO terms related to cell cycle progression and chromatin organization (e.g., GO:0000786, GO:0046982, GO:0007067, GO:0051301). Vascular cell types showed enrichment in hormone signaling (GO:0009734), solute transport (GO:0006820, GO:0008509, GO:0080161), and secondary cell wall biogenesis (GO:0045492). Trichome cells were characterized by terms related to responses to biotic and abiotic stresses, including oxidative stress and defense responses (GO:0031640, GO:0006970, GO:0006805). Guard cells exhibited enrichment in fatty acid biosynthesis (GO:0006633), carboxylic ester hydrolase activity (GO:0052689), and lipid metabolic processes (GO:0044255). Notably, epidermal pavement cells displayed enrichment in photosynthesis-associated categories, including “thylakoid membrane” (GO:0042651), “chloroplast thylakoid membrane” (GO:0009535), “photosynthesis, light harvesting” (GO:0009765), and “response to far-red light” (GO:0010218), in addition to cell wall organization (GO:0009505). This strongly suggests that presence of transcriptionally active nuclear-encoded genes involved in chloroplast function in epidermal pavement cells of R. sericea.

hdWGCNA was applied to systematically uncover coordinated gene expression programs across all cell types. Using 7,574 DEGs, we identified eight distinct co-expression modules (R-se1 to R-se8) characterized by unique expression profiles (Figures 3A–D). Among these, the gray module comprised unassigned genes, whereas the turquoise module was the largest (1,687 genes) and the pink module was the smallest (53 genes). Correlation analysis between modules and cell-type abundance revealed specific associations. For instance, the yellow, brown, green, and pink modules were strongly enriched in meristematic cells, whereas black and blue modules correlated with guard cells, and the red module was specific to xylem cells. Interestingly, the pink module was simultaneously associated with both meristematic and proliferating cells, suggesting overlapping transcriptional programs regulating early growth. In contrast, the turquoise module exhibited no strong correlation with any single cell type (Figure 3E), which may indicate that these are housekeeping or universally regulated genes. These results suggest the presence of both cell-type-specific and cross-lineage regulatory modules. Modules enriched in proliferative or stress-related cell types may represent key transcriptional programs driving adaptive leaf development in R. sericea.

To further elucidate regulatory processes at the species level, we performed GSEA on hub genes identified from the hdWGCNA modules. A total of 724 hub genes were analyzed (Supplementary Table S4). The analysis revealed a significant enrichment of the SPL7_TARGET_GENES set (p-value = 1.548e-08), with an enrichment score (ES) of 0.45 and a normalized enrichment score (NES) of 2.94. This set contained 64 genes, including 57 core-enriched members that contributed strongly to the signal (Figures 4A, B). Squamosa promoter-binding protein-like 7 (SPL7), a member of the SBP transcription factor family, is a pivotal regulator of copper homeostasis in plants (Yan et al., 2017; Mermod et al., 2019).

To reconstruct the developmental trajectory of leaf cells, we performed pseudotime analysis across all captured cell types by Monocle 2. Meristematic cells were designated as the root of the trajectory, representing the least differentiated state. Differentiation was visualized along a continuum from initiation (blue zone) to completion (red zone) (Figure 5A).

Our analysis revealed two distinct regions of differentiation initiation: one originating from the meristematic and proliferating cell clusters, and the other from the epidermal cell population. These dual initiation zones suggest that both internal tissues and epidermal lineages harbor high developmental potential. A striking finding was that the initiation zone of differentiation within the epidermis corresponded to the guard cell cluster (Figure 5B). This observation may suggests the presence of a lineage within guard cells sharing characteristics with protodermal or guard mother cells, exhibiting high developmental potential. Moreover, mesophyll cells were predominantly localized at the terminal end of the trajectory, consistent with their status as fully differentiated, photosynthetically specialized cells.

To investigate the developmental dynamics of trichome lineages, we conducted pseudotime analysis of epidermal cells, including pavement, guard, and trichome cells. we designated meristematic cells as the initial differentiation state and reconstructed the epidermal developmental trajectory. The results showed that the pseudotime trajectory contained a single branch point that separated the epidermal population into three distinct states (Figures 5C, D). As expected, meristematic cells population dominated early differentiation and were enriched in the pre-branch state (state 2). After bifurcation, guard cells and epidermal pavement cells diverged into two terminal fates, occupying states 1 and 3, respectively. Interestingly, trichome cells shared a differentiation trajectory with epidermal pavement cells, but were only temporally later than the latter (Figures 5D, E).

Branch-dependent analysis revealed that the top 100 DEGs across the three states could be grouped into six expression clusters (Figure 5F). Genes in cluster 2 were predominantly associated with trichome differentiation and enriched in defense-related processes, such as “killing of cells of other organisms” (GO:0031640) and “chitin catabolic process” (GO:0006032) (Figure 5G). These results suggest that trichome commitment is accompanied by the activation of defense-associated transcriptional programs, reflecting the dual role of trichomes in morphogenesis and stress adaptation (Viterbo et al., 2002).

Among the top 100 genes showing the most pronounced transcriptional changes along the epidermal differentiation trajectory, we focused on RO6G37307 (a 495-bp gene encoding a 164-amino-acid protein, Supplementary Table S5) as a candidate regulator because it showed markedly specific and elevated expression in the trichome cells cluster (Figure 5H, p-value = 7.6e-161). Furthermore, during trichome differentiation, its expression pattern was highly coordinated with gene sets associated with trichome biological functions identified based on homologous genes(Figure 5F), and its expression timing closely aligned with the differentiation trajectory of trichome cells (Figure 5I), suggesting its potential involvement in the trichome development process.

To explore its potential biological function, we ectopically overexpressed RO6G37307 in A. thaliana. Overexpression markedly increased trichome density on leaf surfaces—approximately sevenfold higher than that of Col-0 wild type (71 vs. 10 trichomes per unit area)—and also altered trichome morphology, with a subset exhibiting four-branched structures rather than the typical three-branched form (Figures 6A–E).

Phylogenetic analysis based on homologous amino acid sequences placed RO6G37307 within a well-supported clade (bootstrap = 98%) alongside Phylloplanin homologs from Rosa and Fragaria, supporting its membership in the Phylloplanin family (Figure 6F). Signal peptide prediction indicated a high-confidence N-terminal cleavage site between residues 27 and 28 (probability = 0.98; Figure 6H). Transmembrane topology analysis suggested a “Globular + SP” architecture, lacking transmembrane helices and carrying a cleavable signal peptide typical of secreted, soluble proteins (Figure 6G). Finally, subcellular localization assays through GFP fusion experiments showed that RO6G37307 protein is targeted to the extracellular space, and its fluorescence signal is found on the cell wall (Figure 6I), consistent with its classification as a secreted phylloplanin-like protein. Together, these results support the role of RO6G37307 as a novel regulator of trichome morphogenesis in R. sericea and point to its potential contribution to epidermal defense adaptation.

Elucidating plant developmental programs is essential for understanding evolutionary processes and adaptive strategies (Li et al., 2024). SnRNA-seq provides unprecedented resolution for analyzing cellular heterogeneity during growth and development (Zhu et al., 2022; Ye et al., 2023; Conte et al., 2024), and has revolutionized developmental biology (Nakayama et al., 2022). Applying snRNA-seq to autumn leaves of R. sericea offers a unique opportunity to reveal the cellular and molecular mechanisms underlying plant responses to plateau environmental warming.

In this study, we leveraged prior knowledge of Arabidopsis marker genes and required that each cluster be annotated with multiple validated marker genes whenever possible, in order to reduce the risk of misannotation through cross-validation. Additionally, we integrated data on the biological functional profiles of the cell clusters and types, and cross-species transcriptomic comparison to provide complementary support for cell type annotation. Using this strategy, we provisionally identified and annotated 11 distinct cell types in R. sericea leaves, resulting in a draft cell atlas encompassing major tissues, including leaf epidermis, mesophyll, and vascular systems. The consistency of our annotations with previous single-cell studies further supports their biological plausibility and demonstrates the feasibility of scRNA-seq in non-model alpine plants (Wang et al., 2022a; Zang et al., 2023; Guo et al., 2025). Importantly, this work fills a gap by extending single-cell technologies to autumn shoots of R. sericea, thereby providing a valuable resource for exploring adaptive mechanisms in alpine ecosystems. However, Cell type annotation in non-model plant species remains challenging due to the scarcity of species-specific markers and high-quality reference atlases. In this study, cell identities in R. sericea were inferred primarily through conserved orthologous markers from A. thaliana, complemented by functional enrichment analyses. Although this integrative strategy reduces misannotation risk through cross-validation, it is inherently indirect and inference-based. In addition, snRNA-seq may underrepresent cytoplasmic transcripts, potentially limiting the detection of certain cell-type-specific expression features, while cross-species marker transfer introduces uncertainty due to evolutionary divergence in gene regulation. Therefore, cross-species correlation analysis should be viewed as an auxiliary consistency check rather than a quantitative measure of annotation accuracy. Accordingly, the cell type assignments presented here should be regarded as putative rather than definitive. Future studies combining reference-based cross-species mapping approaches with independent spatial validation methods, such as RNA fluorescence in situ hybridization (FISH), will be important for further improving annotation confidence and resolution in R. sericea and other non-model plant systems.

Each plant cell type is not only distinct in identity but also exhibits functional heterogeneity and intricate interconnectedness (Mathe et al., 2021; Thibivilliers and Libault, 2021; Amini et al., 2023). These attributes are fundamental to a plant’s ability to adapt to diverse environmental challenges, sustain normal physiological and developmental processes, and coordinate complex biological activities. Examining the functional dynamics of cell populations under specific environmental conditions provides critical insights into adaptive strategies (Mackenzie and Mullineaux, 2022; Roman et al., 2022; Yu et al., 2024).

R. sericea initiates autumnal shoot growth at favorable times, thereby undergoing a secondary phase of vegetative development that effectively extends the growing season. Functional enrichment analyses revealed that autumn leaves of R. sericea retained diverse and indispensable biological functions (such as cell cycle progression, photosynthesis, hormone signaling, metabolic processes and defense responses) necessary for maintaining physiological homeostasis, while also exhibiting robust differentiation dynamics (Supplementary Figure S4). Previous studies suggested that delayed senescence increases root biomass, though not branch biomass (Zohner et al., 2021). However, it remains uncertain whether extended autumnal growth contributes to greater assimilate accumulation in this species. We speculate that, as a deciduous shrub, R. sericea may allocate assimilates gained from autumnal shoots into root storage, potentially forming a reservoir that supports germination and growth in the following spring.

Remarkably, epidermal pavement cells of R. sericea harbor transcriptionally active nuclear-encoded genes involved in chloroplast function. Prior studies have demonstrated that chloroplasts in pavement cells not only perform photosynthesis (Barton et al., 2016), but also contribute to immune responses (Irieda and Takano, 2021; Irieda, 2022). Thus, we propose that epidermal pavement cells may represent an adaptive strategy to optimize both light utilization and defense in the high-altitude environment. By extending light capture and primary photochemistry to the leaf surface, epidermal pavement cells can supplement mesophyll carbon gain during periods of low temperature and limited stomatal conductance, when mesophyll photosynthesis is constrained. However, we still lack direct evidence to prove whether epidermal pavement cells of R. sericea contain functional chloroplasts. Visualization is recommended in future work, such as chlorophyll autofluorescence, confocal imaging and subcellular localization.

While traditional gene expression studies often focus on individual genes, overlooking their coordinated interactions, we employed co-expression network analysis by hdWGCNA to delineate regulatory relationships (Langfelder and Horvath, 2008; Cao et al., 2023). Cell type–phenotype associations across eight defined modules revealed robust correlations between specific clusters and transcriptional programs, and cleanly onto specific biological processes. These functional hubs may reflect key transcriptional landscape that are associated with leaf developmental patterns potentially relevant to alpine adaptation in R. sericea. Collectively, these modules suggest coordinated biological programs rather than isolated gene effects. To link modules to ecological outcomes, targeted validation is needed in the future, such as cell-specific perturbations (RNAi/CRISPR), spatially resolved transcriptomics/metabolomics, and physiological assays.

At the gene-set scale, we further identified physiological adaptation mechanisms. Notably, enrichment of SPL7 target genes—SPL7 being a central regulator of copper deficiency responses—indicates that copper homeostasis plays a potential role in autumn leaf growth of R. sericea under copper-deficient soil conditions in its habitat (Zhang et al., 2002; Bing et al., 2016). Mechanistically, SPL7 enhances copper uptake and utilization by activating copper-responsive genes under deficiency, thereby maintaining micronutrient equilibrium. Our data suggest that the activation of the SPL7-mediated copper homeostasis network is associated with the adaptive response in R. sericea, and may represents a promising candidate mechanism. Although this study suggests a potential key role for the SPL7 module in the leaf development of R. sericea, its direct physiological function and causal relationship require further confirmation through future genetic experiments combined with physiological phenotype measurements.

Analyzing leaf developmental dynamics allows us to elucidate how plants regulate morphogenetic programs to generate specialized functional structures. Such analyses not only deepen our understanding of fundamental developmental mechanisms but also reveal how plants adapt to environmental stress.

Previous pseudotime studies suggested that leaf primordia represent the main point of origin, although the timing of terminal differentiation varies among species (Bai et al., 2022; Wang et al., 2022a; Zang et al., 2023; Liang et al., 2025). Interestingly, work on peanut leaf morphogenesis revealed that epidermal cell fate can be specified even prior to primordium formation (Liu et al., 2021). In our study, the developmental trajectory began with clusters of protoblasts exhibiting strong proliferative and differentiative potential—namely, meristematic and proliferating cells—and culminated in the differentiation of mesophyll cells (Figure 5A). We infer that a hierarchical program in which differentiation potential is progressively restricted and metabolic specialization is sequentially acquired. The positioning of mesophyll cells at the terminal state may imply that their differentiation requires extended developmental input, particularly the coordination of chloroplast biogenesis, cell expansion, and establishment of photosynthetic machinery. This temporal delay likely ensures that mesophyll maturation is synchronized with leaf expansion, thereby optimizing carbon assimilation once the photosynthetic surface has fully developed. Together, these findings describe the diversity of developmental programs shaping leaf formation and suggest that R. sericea may exhibit a trajectory optimized for energy acquisition in the alpine environment, pending further validation.

As key adaptive structures on leaf surfaces, epidermal trichomes serve as a valuable model system for studying plant development, with significance in both fundamental research and commercial applications (Dong et al., 2023), and are evolutionarily conserved across diverse plateau plant taxa (Wu et al., 2024a). During leaf development of R. sericea, epidermal trichomes contribute substantially to resistance against biotic and abiotic stresses. Their presence is associated with the species’ adaptability to high-altitude conditions, reflected by function analyses in cell fate differentiation programs. Based on the temporal regulation of development, trichome cells are closely related to epidermal pavement cells at the late stage of leaf epidermal differentiation. This discovery not only reveals cell fate transitions during epidermal development, thus expanding our understanding of leaf epidermal differentiation, but also suggests that the spatiotemporal initiation of trichome morphogenesis is closely related to environmental adaptation.

Concurrently, our analysis further identified that RO6G37307 is a gene within the Phylloplanin family. Phylloplanins are cysteine-rich, secreted proteins first characterized in Nicotiana tabacum (Shepherd et al., 2005), where they mediate chemical defense on phylloplane surfaces by inhibiting airborne pathogens (Kroumova et al., 2007, 2013; Sahoo et al., 2014; Freire et al., 2017). Beyond this canonical role, our findings suggest the potential functional scope of phylloplanin-like gene by linking the tissue-specific expression of RO6G37307 to trichome morphogenesis. It is important to note, however, that while phylogenetic, structural prediction, and subcellular localization analyses converge to support its identification as a secreted phylloplanin-like protein, they do not constitute definitive proof. A conclusive functional assignment will therefore require targeted biochemical or genetic validation.

In addition, TCP4 and TTG2, known as trichome-specific markers (p-value = 4.4e-29, p-value = 3.6e-19), exhibit unique regulatory timing and high expression levels in R. sericea trichome clusters (Supplementary Figure S8-S11). Previous studies have shown that TCP4 acts as a negative regulator of trichome development: tcp4 mutants display increased trichome density and four-branched morphology (Vadde et al., 2018; Wang et al., 2022c). TTG2 expression in late-stage trichome morphogenesis is co-regulated by the TTG1-GL1-GL3 complex and (mitogen-activated protein kinase) MAPK signaling (Pattanaik et al., 2014; Ma et al., 2024). Upon induction, TTG2 directly activates cytoskeleton-associated genes (e.g., BRICK1), remodeling microtubule and actin networks to promote morphogenesis (Liu et al., 2024). These established roles align with the trichome phenotypes observed upon ectopic overexpression of RO6G37307. Integrative analysis leads us to propose one speculative model wherein RO6G37307 may be associated with leaf epidermal trichome development via interactions with the TTG2 and TCP4 regulatory networks (Figure 7). In this model, RO6G37307 might affect extracellular ROS homeostasis, potentially influencing the MAPK cascade to modulate TTG2 expression and TCP4 activity (Xu et al., 2013; Taj et al., 2014; Rayapuram et al., 2018), thereby creating a permissive environment for trichome morphogenesis (Ma et al., 2024). However, this proposed pathway is inferred from transcriptomic correlations and published literature rather than demonstrated causal relationships. Accordingly, further experimental evidence will be required to test and refine this hypothetical regulatory framework. At present, our findings suggest a potential connection between RO6G37307 expression, leaf epidermal trichome cell differentiation, and stress-related signaling processes, which may contribute to adaptive traits in R. sericea under alpine environments, but such adaptive significance remains to be directly demonstrated.