Single-cell transcriptomics reveals cellular and genetic mechanisms of alpine adaptation in Rosa sericea
Hengning Deng, Jian Ru, Zhenlong Liang, Zhongyu Tang, Yang Wang, Wenqin Yuan, Liangying Li, Yu Feng, Xinfen Gao
TL;DR
This study uses single-cell RNA sequencing to uncover how Rosa sericea adapts to high-altitude environments at the cellular and genetic levels.
Contribution
The first single-cell transcriptional atlas of Rosa sericea leaves, revealing adaptive mechanisms to alpine conditions.
Findings
A developmental continuum from cell proliferation to photosynthetic maturation was identified in R. sericea leaves.
Transcriptionally active genes in epidermal cells suggest roles in chloroplast function and copper homeostasis for high-altitude adaptation.
A regulatory module involving RO6G37307, TTG2, and TCP4 is linked to trichome development in alpine adaptation.
Abstract
Plant development is shaped by environmental conditions, and its adaptation to climate change is crucial for biodiversity conservation. The extreme climate of the Qinghai-Tibet Plateau makes it an ideal system for studying plant adaptive strategies. Rosa sericea, a dominant alpine shrub, exhibits remarkable morphological plasticity, but its molecular and cellular adaptation mechanisms are still unclear. In this study, we integrated single-nucleus RNA sequencing (snRNA-seq) with high-dimensional weighted gene co-expression network analysis (hdWGCNA), gene ontology (GO) enrichment, gene set enrichment analysis (GSEA), pseudotime trajectory inference, and gene overexpression techniques to profile 31,796 cells from R. sericea leaves. We constructed a draft single-cell transcriptional atlas with putative annotation of 11 leaf cell types and identified eight co-expression gene modules linked…
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| Cell types | Marker genes | References |
|---|---|---|
| Mesophyll cell |
| ( |
| Epidermis cell |
| ( |
| Bundle sheath cell |
| ( |
| Phloem parenchyma cell |
| ( |
| Phloem cell |
| ( |
| Xylem cell |
| ( |
| Companion cell |
| ( |
| Meristematic cell |
| ( |
| Proliferating cell |
| ( |
| Guard cell |
| ( |
| Epidermal pavement cell |
| ( |
| Trichome cell |
| ( |
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
TopicsPlant Molecular Biology Research · Photosynthetic Processes and Mechanisms · Plant Gene Expression Analysis
Introduction
1
Global climate warming presents significant challenges to plant survival (Kumar et al., 2024a; Seth and Sebastian, 2024; Sun et al., 2025), profoundly affecting their growth and development (Li et al., 2019; Saladin et al., 2020; Reich et al., 2022; Wang et al., 2022b; Sigdel et al., 2024). These effects are especially evident in evolutionary and developmental characteristics (Wright et al., 2017; Yuan et al., 2019; Cabon et al., 2022; Gao et al., 2022; Meng et al., 2024). One of the most notable consequences is the advancement of leaf unfolding in spring and delayed senescence during autumn (Fu et al., 2015; Dow et al., 2022; Marques et al., 2023; Sun et al., 2024). While numerous studies have explored plant adaptation strategies to climate change, most focus on large-scale patterns, such as vegetation (Saladin et al., 2020; Antão et al., 2022; Cabon et al., 2022; Mirabel et al., 2023; Pu et al., 2025) and forest phenology (Piao et al., 2019; Shen et al., 2022), with few investigating species-specific adaptability (Weih, 2009; Repo et al., 2021; Martinez Del Castillo et al., 2022; Zlobin et al., 2024). In-depth research at the species level is essential for understanding how plants adjust their survival strategies to rapidly changing climatic conditions, which is crucial for elucidating their survival mechanisms and ensuring biodiversity conservation (Ram, 2024). As a globally recognized biodiversity hotspot, the Qinghai-Tibet Plateau exhibits heightened sensitivity to climate change (Meng et al., 2023; Tang et al., 2023). Ecological adaptation studies in this region not only enrich plateau ecological theories but also provide a scientific foundation for biodiversity conservation in high-altitude ecosystems (Pu et al., 2025; Sun et al., 2025). Key adaptive strategies in alpine plants include specialized leaf morphology and physiology, such as smaller leaf surface areas, thicker cuticles, and modified stomatal patterns, which reduce water loss and mitigate solar radiation (Jin et al., 2024). The development of glandular and non-glandular trichomes on leaves also plays a vital role in enhancing survival, as these structures can enhance the plant’s ability to withstand various stresses (Dong et al., 2023). Additionally, the trade-off between trichome density and chemical defenses may have significant implications for plant survival strategies in resource-limited environments (Xu et al., 2025). These strategies not only ensure plants survival in harsh alpine conditions but also contribute to biodiversity across the diverse habitats of the plateau. Investigating these adaptations provides insights into how plants face climate change challenges and supports high-altitude ecosystem conservation.
As a prominent shrub of the Qinghai-Tibet Plateau (Gao et al., 2015), Rosa sericea is widely recognized for its extraordinary phenotypic plasticity (Ullah et al., 2022a; Jiao et al., 2023), a key factor in maintaining the ecological balance of the region. During the spring and summer months, R. sericea undergoes vigorous growth, followed by a physiological dormancy phase in autumn. This dormancy acts as an adaptive strategy, shielding the plant from imminent environmental extremes. However, climate change on the Qinghai-Tibet Plateau (Li et al., 2019) has lengthened the period of autumn warmth, disrupting the natural growth patterns of this plant (Wang et al., 2022b; Wu et al., 2022; Sigdel et al., 2024). Recent field studies have shown that R. sericea often initiates a second vegetative growth cycle during the warm autumn, sprouting new shoots in the same year. This characteristic makes R. sericea a potentially excellent model for studying how climate change modifies plant adaptation strategies. Although initial research has focused on leaf morphology and metabolomics (Ullah et al., 2022b; Jiao et al., 2023), a deeper understanding of its developmental mechanisms at the cellular and molecular levels remains lacking, limiting our knowledge of how R. sericea adapts to its environment.
The emergence of single-cell RNA sequencing (scRNA-seq) technology has revolutionized the study of biological phenomena at the resolution of individual cells (Qin and Tape, 2024; Lee et al., 2025). In recent years, scRNA-seq has provided significant insights into the development of complex tissues (Jovic et al., 2022; Wang et al., 2023; Conte et al., 2024; Zhu et al., 2025), greatly enhancing our understanding of cellular identity (Liu et al., 2023; Ye et al., 2024; Guo et al., 2025), organ heterogeneity (Zhang et al., 2023; Marchant and Walbot, 2025), gene regulation (Guo et al., 2025; Liang et al., 2025), and cellular metabolic networks (Li et al., 2023; Wu et al., 2024b). In this study, we applied snRNA-seq to R. sericea tender leaves, mapping its transcriptional landscape and uncovering the cellular-genetic mechanisms underlying leaf adaptive development. To harness this potential at single-cell resolution, we defined the core objectives of this study as follows:(1) How is cellular heterogeneity organized in R. sericea leaves, and what are the developmental lineages of major cell types? (2) Which transcriptional programs and cell types show signatures of being shaped by or responding to high-altitude stresses? (3) What is the regulatory architecture contributing key adaptive traits, such as trichome development, and how does it integrate developmental cues with potential environmental signals? Through this work, we aim to (i) contribute to a cell-resolved understanding of leaf adaptive development in an alpine shrub, (ii) explore potential links between molecular networks and ecologically relevant phenotypes, and (iii) offer testable hypothesis model for future functional and evolutionary studies.
Materials and methods
2
Plant material collection and preparation
2.1
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.
Single-cell sequencing library construction process and cell clustering of R. sericea. (A) Habitat and flowering plants of R. sericea. (B) Tender autumn leaves of R. sericea. (C) Overview of the scRNA-seq library construction process for R. sericea leaves (workflow diagram adapted from https://grcf.jhmi.edu/service/10x-single-cell/). (D) UMAP plot showing the clustering of leaf cells, with each dot representing an individual cell. The 22 distinct colors correspond to 22 specific cell clusters. Replicates are indicated as R-se1, R-se2, and R-se3. (E) Pearson correlation analysis between cell clusters. Each square represents the correlation coefficient between two clusters. Darker colors indicate higher correlation (values closer to 1). (F) Heatmap of the top 10 DEGs across 22 cell clusters. The horizontal axis represents cell populations, and the vertical axis denotes the DEGs for each population. Color intensity reflects gene expression levels, with red indicating higher expression and blue indicating lower expression.
Nucleus isolation protocol
2.2
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.
Construction and sequencing of snRNA-seq libraries
2.3
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.
SnRNA-seq data preprocessing
2.4
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.
Identification of highly variable genes, cell clustering, and visualization
2.5
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).
Tissue-specific marker gene identification and leaf cell-type annotation
2.6
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.
Cross-species mapping to Arabidopsis leaf snRNA-seq reference
2.7
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.
High-dimensional weighted gene co-expression network analysis
2.8
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).
Gene ontology enrichment analysis and gene set enrichment analysis
2.9
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.
Pseudotime trajectory analysis
2.10
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.
Construction of overexpression vector and plant transformation
2.11
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.
Selection of high-expression lines and phenotypic evaluation
2.12
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.
Protein functional analysis
2.13
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.
Results
3
Construction of a high-resolution single-cell transcriptomic atlas of R. sericea leaves
3.1
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).
Classification and identity annotation of leaf cells
3.2
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).
Cell-type annotation of R. sericea leaves. (A, B) UMAP embeddings of all single cells; each point represents one cell and is colored by the assigned cell type. (C) Bar chart showing cell counts per type; colors match the palette in panels. (D, E) Dot plots of canonical marker genes across cell types; dot size indicates the fraction of cells expressing each gene, and color encodes the average scaled expression (warmer colors indicate higher expression). (F) Subclustering and re-annotation of the epidermal lineage by 15 marker genes, highlighting pavement, guard, and trichome subtypes on UMAP.
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.
Cross-species comparison supports the overall plausibility of cell type annotation
3.3
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.
Gene ontology analysis reveals functional diversity across cell types
3.4
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.
Gene co-expression network analysis identifies key regulatory modules
3.5
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.
*The hdWGCNA analysis of R. sericea leaf cells. (A) Dendrogram of gene clustering. Each branch represents an individual gene, with color coding at the bottom indicating the gene’s association with specific co-expression modules. The “grey” module represents genes that were not assigned to any specific co-expression module. (B) UMAP projections of the harmonized module eigengenes (hMEs) for the eight modules. (C) Correlation analysis between the eight modules (*p < 0.05, **p < 0.01, **p < 0.001). (D) UMAP plot displaying the top hub gene from each module. (E) Bubble plots showing the performance of hMEs from each module across different cell types. The dot plot illustrates the average expression of module-specific Module Eigengenes (ME) across different cell clusters, ranked by their kME values within each module.
GSEA suggests a potential association between SPL7-related transcriptional programs and copper homeostasis
3.6
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).
GSEA analysis based on target gene sets during R. sericea leaf development. (A) Gene set enrichment plot. The upper section displays the enrichment score. Vertical lines in the middle indicate the position of each gene within the ranked gene set list. The lower section shows the distribution of gene ranks, represented as a gray area plot. (B) The top 3 gene sets with the highest enrichment scores. Gene set names, p-values, and adjusted p-values are shown in the upper right corner.
Global pseudotime analysis of leaf development
3.7
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).
Pseudotemporal differentiation and development of R. sericea leaves. (A) Global pseudotemporal annotation of leaf cells, mapping the “pseudo-time” of Monocle 2 onto the two-dimensional UMAP space. (B) Pseudotemporal annotation of leaf epidermal cells, mapping the “pseudo-time” of Monocle 2 onto the two-dimensional UMAP space. (C) Ordering of leaf epidermal cells along pseudotemporal differentiation trajectories. (D) States of leaf epidermal cells along a pseudotemporal differentiation trajectory. (E) Distribution of four cell types in the leaf along pseudotemporal differentiation trajectories. (F) Heatmap showing the top 100 genes with the most significant changes along the differentiation trajectory. The branch point in the middle indicates the start of pseudo-time. The topology used to visualize pseudotemporal development in the four cell types is outlined with red dashed line in the figure. Key genes involved in trichome differentiation are shown on the right, with genes highlighted in bold being the focus of this study. (G) GO functional enrichment of trichome fate-determining genes among the top 100 genes along the leaf epidermal differentiation trajectory. (H) The violin plot of RO6G37307 gene expression in R. sericea leaf epidermal cell subtypes. (I) Expression trends of RO6G37307 in relation to cell fate determination. The dotted line indicates the expression level changes of RO6G37307 along the pseudo-temporal differentiation trajectory of trichomes.
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.
Pseudotime trajectory analysis of trichome differentiation
3.8
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).
RO6G37307 positively regulates trichome development in leaf epidermis
3.9
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).
*Functional validation of the RO6G37307 gene. (A) Relative expression levels of RO6G37307 in overexpression (OE) lines. AtActin was used as the internal reference gene for RT-qPCR analysis. Data represent the mean ± SD from three biological replicates, and the relative expression levels are shown above the bars (***p < 0.001). (B) Comparison of trichome density on leaf surfaces between the OE line RO6G37307–2 and wild-type (Wt). Data represent the mean ± SD from three biological replicates, and the gray dots represent the recorded data points. The inset shows the three counting areas per leaf (boxed regions) (**p < 0.001). (C) Trichome phenotypes on rosette leaves of the OE line RO6G37307-2. Red arrows indicate four-branched trichomes. (D) Representative trichome types on adaxial leaf surfaces of OE line RO6G37307–2 under scanning electron microscope. (E) Comparative visualization of adaxial leaf trichomes under a digital camera: Wt (left) and RO6G37307-2 (right). (F) Phylogenetic analysis of RO6G37307 homologs. Bootstrap support values >50% (from 5,000 replicates) are shown on the branches of the Neighbor-Joining tree. (G) Transmembrane domain prediction for RO6G37307. Top: Amino acid sequence of transmembrane helices (orange) and extracellular regions (blue). Bottom: Confidence scores. (H) Signal peptide prediction. The blue dotted line indicates the cleavage site at the N-terminus. (I) Subcellular localization prediction. In the control, free GFP displayed a uniform distribution in the cell periphery and nucleus. In contrast, RO6G37307-GFP localized to the cell periphery prior to plasmolysis. After plasmolysis, the GFP signal remained at the cell wall (gree) rather than retracting with the plasma membrane (red, as indicated by the arrow), supporting its extracellular secretion and cell wall localization.
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.
Discussion
4
SnRNA-seq establishes a cellular developmental atlas of a plateau plant in autumn
4.1
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.
Functional heterogeneity and physiological homeostasis of leaf cells underpin adaptability in R. sericea
4.2
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.
Developmental trajectories provide insights into plateau plant adaptation
4.3
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.
RO6G37307 regulates morphogenesis and functional acquisition of leaf epidermal trichomes
4.4
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.
A proposed mechanism for R. sericea leaf trichome differentiation. RO6G37307-mediated regulation of epidermal trichome development and functional acquisition (↑/↓: up/down-regulation of expression).
Conclusion
5
This study used single-cell transcriptomics to construct the cell atlas and characterize the temporal program of autumn leaves development in R. sericea. Our analyses highlight several features that may be relevant to alpine adaptation, including presence of transcriptionally active nuclear-encoded genes involved in chloroplast function in epidermal pavement cells, transcriptomic signatures consistent with SPL7-linked copper homeostasis, and a putative RO6G37307–TTG2–TCP4 regulatory module associated with trichome development. These strategies could collectively contribute to the plant’s resilience to plateau stress. Overall, our findings offer insights into potential molecular mechanisms, involving developmental and physiological adjustments, that may contribute to the adaptation of R. sericea to high-altitude environments. These observations advance our understanding of alpine plant adaptation and establish a foundation for future research into adaptive strategies in other species, with potential implications for efforts to improve plant resilience under climate change.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Amini S. Doyle J. J. Libault M. (2023). The evolving definition of plant cell type. Front. Plant Sci. 14. doi: 10.3389/fpls.2023.1271070, PMID: 37692436 PMC 10485272 · doi ↗ · pubmed ↗
- 2Antão L. H. Weigel B. Strona G. Hällfors M. Kaarlejärvi E. Dallas T. . (2022). Climate change reshuffles northern species within their niches. Nat. Climate Change 12, 587–592. doi: 10.1038/s 41558-022-01381-x · doi ↗
- 3Bai Y. Liu H. Lyu H. Su L. Xiong J. Cheng Z. M. (2022). Development of a single-cell atlas for woodland strawberry (Fragaria vesca) leaves during early Botrytis cinerea infection using single cell RNA-seq. Hortic. Res. 9, uhab 055. doi: 10.1093/hr/uhab 055, PMID: 35043166 PMC 8969069 · doi ↗ · pubmed ↗
- 4Barton K. A. Schattat M. H. Jakob T. Hause G. Wilhelm C. Mckenna J. F. . (2016). Epidermal pavement cells of Arabidopsis have chloroplasts. Plant Physiol. 171, 723–726. doi: 10.1104/pp.16.00608, PMID: 27288524 PMC 4902630 · doi ↗ · pubmed ↗
- 5Becht E. Mc Innes L. Healy J. Dutertre C. A. Kwok I. W. H. Ng L. G. . (2018). Dimensionality reduction for visualizing single-cell data using UMAP. Nat. Biotechnol. 37, 38–44. doi: 10.1038/nbt.4314, PMID: 30531897 · doi ↗ · pubmed ↗
- 6Bing H. Wu Y. Zhou J. Li R. Luo J. Yu D. (2016). Vegetation and cold trapping modulating elevation-dependent distribution of trace metals in soils of a high mountain in eastern Tibetan Plateau. Sci. Rep. 6, 24081. doi: 10.1038/srep 24081, PMID: 27052807 PMC 4823730 · doi ↗ · pubmed ↗
- 7Cabon A. Kannenberg S. A. Arain A. Babst F. Baldocchi D. Belmecheri S. . (2022). Cross-biome synthesis of source versus sink limits to tree growth. Science 376, 758–761. doi: 10.1126/science.abm 4875, PMID: 35549405 · doi ↗ · pubmed ↗
- 8Camacho C. Coulouris G. Avagyan V. Ma N. Papadopoulos J. Bealer K. . (2009). BLAST+: architecture and applications. BMC Bioinf. 10, 421. doi: 10.1186/1471-2105-10-421, PMID: 20003500 PMC 2803857 · doi ↗ · pubmed ↗