In June 2018, samples of R. sanguineum var. haemaleum for transcriptome sequencing were collected between 11:00am and 1:30pm within two successive days from four sites along an elevational gradient ( Figure 1 ) in the northern Gaoligong Mountains (N 27°49’53”, E 98°27’03” - N 28°05’13”, E 98°45’31”) which is located in southern of Hengduan Mountains (Liu et al., 2022a). Three individuals were sampled across a small range of 20 m² from each site. The entire late flower bud and entire early leaf bud were sampled from each individual, immediately frozen in liquid nitrogen in the wild, and preserved in a -80°C ultra-low temperature freezer prior to processing for total RNA isolation. Vouchers of each individual plant sampled were collected and deposited in the Herbarium of Kunming Institute of Botany (KUN), Chinese Academy of Sciences.

Total RNA was extracted and purified separately from each organ using a Spectrum TM Plant Total RNA Kit (STRN250, Sigma) according to the manufacturer’s protocols. The RNA purity and quality were assessed with a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was evaluated through agarose gel electrophoresis. Three biological replicates (from three plants) of each flower bud and leaf bud material were included for each elevational site. A single cDNA library was separately constructed for each flower bud and each leaf bud of each individual from each site, and subsequently sequenced on an Illumina HiSeq X Ten sequencing platform (San Diego, CA, USA), generating approximate 6 Gb paired-end reads (2 × 150 bp) for each sample. Both cDNA library preparation and Illumina sequencing were carried out by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China).

We first cleaned the raw sequence reads by using Trimmomatic v. 0.38 (Bolger et al., 2014) with default settings to remove low-quality reads. We further evaluated the quality of the remaining reads using FastQC (Andrews, 2010). All subsequent analyses were based on these cleaned reads. A reference transcriptome was generated by de novo assembly of the combined clean reads (leaf buds, flower buds, in three replicates from each gradient) performed using Trinity v. 2.6.5 (Haas et al., 2013) with default parameters. Assembly statistics were obtained using the TrinityStats.pl script in the Trinity package. The longest transcripts were considered as the non-redundant unigenes. We also used HISAT2 v. 2.1.0 (Kim et al., 2019) to assess assembly quality, by mapping reads back to the assembled transcripts to count the overall alignment rates.

To lower the redundancy in the dataset, low-coverage artifacts or redundancies were removed using the CD-HIT-EST v. 4.7.0 (Fu et al., 2012) with setting word length to 10 and an identity threshold of 90%. Downstream analyses were performed on the final filtered transcripts. To determine the transcriptome completeness of the assembly, Benchmarking Universal Single-Copy Orthologs tools (BUSCO, v. 4.0.6) (Simao et al., 2015) were used to obtain the percentage of single-copy orthologs represented in the embryophyte database (odb10, 1,614 single-copy orthologues) and also to evaluate the completeness of transcript assemblies. To annotate the assembled unigenes, we downloaded a protein reference of R. delavayi, a species closely related to R. sanguineum from the whole genome sequencing project deposited in GigaDB (Zhang et al., 2017). Open reading frames (ORFs) were first predicted from each filtered assembled transcript using TransDecoder v. 5.5.0 (Haas et al., 2013) with a minimum length of 100 amino acids, and the predicted ORFs were scanned to find homology profiles with a cut-off e-value of 1e-10 against the reference protein database. Unigenes and predicted protein sequences were used as queries to search protein databases using the BLASTP v. 2.5.0, setting the e-value cutoff to 1e-10. Queries were performed against the NCBI non-redundant (NR) and UniProtKB/Swiss-Prot databases. We also performed additional functional annotations with DIAMOND (Buchfink et al., 2015) hits against the eggNOG database (Huerta-Cepas et al., 2017), which summarized available functional information from the different proteins databases, including GO (Gene Ontology), COGs/KOGs (Clusters of Orthologous Groups, containing both prokaryotic and eukaryotic clusters), Pfam (Protein families) and KEGG (Kyoto Encyclopedia of Genes and Genomes) (Hernández-Plaza et al., 2023). The best hit was used as the final annotation.

To quantify transcript abundance, we applied the alignment-based methods by mapping all of the cleaned reads from the flower bud samples separately for each biological replicate back to the non-redundant unigenes of the assembled reference transcripts using RSEM v. 1.3.1 (Li and Dewey, 2011) and Bowtie2 v. 2.3.5 (Langmead and Salzberg, 2012) for alignment. When the transcript abundance for each biological replicate flower bud sample had been obtained, we generated a gene expression matrix that was constructed from a matrix of read count with a Trinity script. The differentially expressed unigenes (DEUs) analysis were performed with the DESeq2 package (Love et al., 2014) among the four elevation gradients. Samples from the three high elevations were separately compared to those harvested at lower elevation (3,000 m) to construct comparison groups. DEUs were considered those with false discovery rate (FDR) adjusted p values ≤ 0.05 and absolute values of log2 (fold change) ≥ 1. To compare gene expression values across four populations at different elevations, we used the trimmed mean of M-values normalization (TMM), as implemented in the R package edgeR. All downstream analyses were implemented based on the normalized expression data matrix (TMM normalization). The functional enrichment analyses of DEUs from each comparison were further processed as described previously (Ye et al., 2021), including GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes). GO terms and metabolic pathways with p values ≤ 0.05 were considered significantly enriched by DEUs.

WGCNA (Weighted gene co-expression network analysis) is an analytical method to explore the relationship between modules and concerned phenotypes and to mine key genes in co-expression networks (Langfelder and Horvath, 2008). To further identify the relevant regulatory gene modules related to elevation, a gene co-expression network was constructed on flower and leaf bud organs, respectively, using the WGCNA package in R (Langfelder and Horvath, 2008). This pipeline included principally three steps. Firstly, a hierarchical cluster analysis was carried out after removing low correlation genes, or those with low expression levels using the hclust function. Secondly, to meet the prerequisite of scale-free network distribution, the soft-power threshold β was determined by the function “sft$powerEstimate” to create an adjacency matrix. To better assess the correlation-based association among gene expression patterns, the adjacency matrix was further converted to a topological overlap matrix (TOM), and a gene connection network was constructed. Finally, gene modules were identified and clustered by a dynamic tree cut method based on the eigengenes (ME) of each module, and the modules with closer distances were merged to obtain the appropriate modules. For the two different organs, the threshold to merge similar modules (mergeCutHeight) and minimal gene module size (minModuleSize) were all set to 0.2 and 100, respectively, except the expression split threshold (deepSplit, flower bud was 2 and leaf bud was 1).

To further investigate the gene modules associated with different environmental conditions, the correlation coefficients between module eigengenes and different altitudes were calculated and statistically tested using Pearson correlation analyses. Modules with absolute values of correlation coefficients > 0.60 (P < 0.01) were considered as altitude-relevant modules: the larger the absolute values of correlation coefficient, the higher the correlation between the module and altitude. Functional enrichment analyses of DEUs in altitude-related modules were carried out by using the AnnotationForge and clusterProfiler package in R (Yu et al., 2012; Carlson and Pages, 2018).

The transcriptome data of 24 cDNA libraries originating from the two organs (late flower buds and early leaf buds) for three individuals from each of the four elevational samples were obtained in this study. After trimming adapter sequences and removing low-quality sequences, we generated 40,622,598 – 59,403,462 sequence reads with the Q30 percentage of 89.22% – 93.82% and the GC percentage of 46.35% – 48.84% ( Supplementary Table S1 ). A total of ~174.4 Gb clean data was generated after merging all reads of the 24 libraries. The Q30 (99.9% base call accuracy), GC-content of the merged data were 92.60% and 43.52%, respectively ( Table 1 ). All of the clean reads were assembled into transcripts using the de novo assembly tool. In total, 335,121 transcripts were identified with the average length of 1,156 bp and N50 length of 1,707 bp. We selected the longest transcript as the representative for each cluster. All the transcripts were further filtered to obtain 162,995 unigenes with the average length of 907 bp, N50 length of 1,201 bp and GC-content of 42.02% ( Table 1 ). The majority of the assembled unigenes were in the ranges of 300–2,000 bp and 14,698 (9.02%) were over 2,000 bp ( Table 1 , Supplementary Figure S1B ). Assessment of the transcriptome completeness identified 1,543 (95.60%) complete and fragmented BUSCOs based on the 1,614 conserved BUSCO embryophyte orthologs ( Supplementary Figure S1A ). The reads of the 24 cDNA libraries were aligned by mapping the reads back to the assembled reference transcript, with mapping rates ranging from 93.85% to 96.07% ( Supplementary Table S1 ). These indicated that our transcriptomes were well assembled, and were of good quality.

To determine putative functions of the identified unigenes of R. sanguineum var. haemaleum, partial and complete open reading frames (ORFs) were predicted using the TransDecoder software, and CD-HIT-EST was used to cluster protein coding transcripts to obtain the non-redundant proteins. Finally, a total of 64,233 coding proteins were translated from the filtered assembly ( Table 1 ). Similarity searches for each unigene against public databases, including National Center for Biotechnology Information non-redundant (NR) protein database, UniProt, COG and Pfam found 49,7769 (77.49%), 29,694 (46.23%), 42,709 (66.49%), 44,380 (69.09%) matches, respectively. Taken together, 53,334 (83.03% of the total) of the unigenes showed a match in at least one of these searchable databases, and 27,802 (43.28%) showed significant matches in all of the four databases ( Table 1 ; Supplementary Figure S1C ). Furthermore, 22,454 (34.96%) unigenes were classified into different GO functional groups. The most abundant GO subcategories for biological processes, molecular functions, and cellular components were cellular process 17,635 (27.45%), catalytic activity 10,098 (15.72%) and cell 19,384 (30.18%), respectively ( Table 1 ; Supplementary Figure S2A ). We assigned 15,553 (24.21%) unigenes to 405 KEGG metabolic pathways ( Table 1 ), which represents a valuable resource for investigating specific processes, functions and pathways. Details regarding the KEGG annotation can be found in Supplementary Figure S2B .

All expression analyses were performed separately for the three biological replicates of the flower and leaf organs. A gene expression matrix with 24 columns and 162,995 lines were generated. Each column represented a sample and each line corresponded to the expression of a unigene. The gene expression density and gene expression distribution of 24 individual libraries in late flower buds and early leaf buds can be found in Supplementary Figure S3 . PCA analysis and correlation matrix showed a good correlation between the replicate sets for each elevation of the 24 individual libraries, each elevation and different organ clustered together, showing a population-specific and organ-specific pattern ( Figure 2A ), which indicated that the RNA-seq libraries were reliable.

Based on the expression profiles at different elevational sites, differentially expressed unigenes (DEUs) were identified by comparing them to the lowest elevation (3,000 m). A total of 2,358 DEUs (up-regulated: 1,084, down-regulated: 1,274) were obtained in the comparison 3,800 m vs. 3,000 m, followed by 3,600 m vs. 3,000 m (2,294, up-regulated: 1,027, down-regulated: 1,267) and 3,200 m vs. 3,000 m (806, up-regulated: 349, down-regulated: 457) in late flower bud, while 2,711 DEUs (up-regulated: 1,629, down-regulated: 1,082), 1,887 DEUs (up-regulated: 1,058, down-regulated: 829), and 821 DEUs (up-regulated: 422, down-regulated: 399) were screened in early leaf bud, respectively ( Figure 2B ). These graphs indicated that the number of DEUs increased with increasing elevation span in each comparison group of the two organs. In total, we generated 130 DEUs in three comparisons of late flower bud ( Figure 2C ) and 102 DEUs in early leaf bud ( Figure 2D ). Details of the DEUs pattern are shown in Supplementary Figure S4 .

The GO and KEGG functional enrichment analyses of the shared DEUs in the three comparisons showed that the GO enrichment terms of the DEUs were classified into three categories of gene ontologies: biological processes (BP), cellular components (CC) and molecular functions (MF). As shown in Supplementary Figure S5A , the most significantly enriched GO terms in the two different organs (late flower bud, early leaf bud) were sesquiterpene metabolic processes (GO:0051761) and negative regulation of leaf senescence (GO:1900056), respectively, while the KEGG enrichment analysis showed that the DEUs were significantly enriched in sesquiterpenoid and triterpenoid biosynthesis (ko00909), fat digestion and absorption (ko04975), photosynthesis-antenna proteins (ko00196) and tropane, piperidine and pyridine alkaloid biosynthesis (ko00960) ( Supplementary Figure S5B ). The genes enriched in these GO terms and KEGG pathways may be closely related to the species’ elevational adaptability.

In the WGCNA analysis, to further explore the specific genes that were highly associated with high-elevation adaptation in different organs, separately performed based on the DEUs from late flower buds and early leaf buds, we independently identified 4,242 and 4,581 DEUs among the three comparisons in the two organs. Based on the scale-free topology criterion (scale-free R ² of 0.80 and mean connectivity close to 0), the soft-power threshold of β =15 and 8 were selected to generate a hierarchical tree of flower and leaf bud, respectively ( Figures 3A, B ). The gene modules were detected based on the topological overlap matrix (TOM), and genes from the two organs were all assigned into 15 distinct modules based on the similarity of their expression patterns. The numbers of the genes in each module varied greatly from 124 to 880 in late flower bud, and 127 to 759 in early leaf bud, regardless of those unclassified genes in grey modules ( Figures 3C, D ).

To characterize the key modules associated with elevation in R. sanguineum var. haemaleum, we calculated the correlation coefficients between the module and elevation. Modules with the module-trait relationships |r| > 0.60 and p < 0.01 were selected as the key ones that were significantly associated with elevation. Among the 15 gene co-expression modules, four specific modules were obtained in different organs of flowers and leaves, and all of them were significantly negatively correlated with altitude (r < 0), implying that DEUs in these modules may be positively expressed at lower elevations in the flower and leaf organs of R. sanguineum var. haemaleum. The turquoise module was significantly correlated with elevation (r = -0.88, p = 2e-04) ( Figure 4 ), followed by the green module (r = -0.86, p = 3e-04), brown module (r = -0.8, p = 0.002) and red module (r = -0.76, p = 0.004), with 880, 134, 386 and 307 DEUs assigned respectively, in late flower bud. While, in early leaf bud, the tan module (r = -0.89, p = 1e-04) with 171 DEUs, the green module (r = -0.87, p = 3e-04) with 397 DEUs, the red module (r = -0.85, p = 5e-04) with 288 DEUs, and the turquoise module (r = -0.76, p = 0.004) with 759 DEUs were significantly correlated ( Figure 4 ).

In the GO and KEGG analyses, conducted to investigate the biological function of the genes in each altitude-related module, the former indicated terpenoid metabolic processes (GO:0006721), response to water deprivation (GO:0009414), plant hormone signal transduction (ko04075) and sesquiterpenoid and triterpenoid biosynthesis (ko00909) as the most significantly enriched functional annotations in late flower buds, and phenylpropanoid biosynthetic processes (GO:0009699), leaf senescence (GO:0010150), flavonoid biosynthesis (ko00941) and photosynthesis-antenna proteins (ko00196) as the most significantly enriched terms/metabolic pathways in early leaf buds. These were the key enrichment objects for the subsequent analysis ( Figure 5 ). Additionally, there were some GO categories or metabolic pathways that were shared in the two organs, while the organ-specific enrichment types were the most striking ones worthy of further analysis, such as sesquiterpene biosynthetic processes (GO:0051762) that occur only in floral organ, and flavonoid biosynthesis (ko00941) unique to leaf organ.

Based on the result of the GO and KEGG enrichment of DEUs, we confirmed the functional annotations of all candidate genes ( Table 2 ). The sample expression patterns were clustered and visualized in a heatmap to clearly understand the expression of all candidate genes from the two organs ( Supplementary Figure S6 ). In both organs, most gene expression patterns were negatively correlated with altitude, in particular genes related to terpene metabolism in the late flower bud (RsHMGR, RsPT, RsTPS) and anthocyanin metabolism (RsCHS, RsF3’5’H, RsOMT) in the early leaf bud ( Table 2 ). The genes in the terpenoid biosynthesis pathway identified in late flower bud and in the anthocyanin metabolism pathway of early leaf bud were found highly expressed at low elevations and decreased with increasing elevation ( Figure 6 ). Additionally, some genes had multiple homologous copies, with the largest number of 8 and 5 in RsTPS and RsCHS respectively ( Table 2 ), which to a certain extent indicated the importance of these two genes in their respective pathways. However, genes related to each node were not enriched across the entire metabolic pathway ( Figure 6 ).

In this study, we carried out a transcriptome expression profiling analysis in flower and leaf bud sampled along an elevational gradient of R. sanguineum var. haemaleum. The results showed a organ-specific expression profile pattern ( Figure 2A ), which has been proven in other studies (Ye and Varner, 1991; Leung et al., 2014; Tu et al., 2021). Our results were similar to those of Hao et al. (2019) whose work based on comparative transcriptomics of three closely related bird species at different elevations. The sample replicates, i.e., sampled from three individual plants at each elevation, clustered together but separately for elevation ( Figure 2A ; Supplementary Figure S6A ), indicating that plants growing in the same habitat have a more similar expression profile, while differing from those growing at different elevations. Environmental conditions vary at different elevations, and factors, such as temperature, moisture, light quality, and oxygen, will likely to be involved in determining gene expression levels (Gibson, 2008; Liu et al., 2022b). In addition, by comparing the high and low elevational samples of different organs, i.e., late flower buds and early leaf buds, the number of significantly DEUs in their respective comparisons showed a positive correlation with the magnitude of the elevational distance, i.e., the greater the elevation span, the greater the number of differentially expressed unigenes ( Figure 2B ). This is consistent with other studies on Potentilla L. (Ma et al., 2015), Bos L. (Ma et al., 2022a) and other species of Rhododendron (Liu et al., 2022b). In this sense, one could hypothesize that along a mountain cline, with increasing elevation, the stresses caused by environmental factors become more severe, and as adaptation of species to these changing environments additional genes are activated, also known as environmental adaptation (Sun et al., 2018; Wang et al., 2021b). Therefore, elevation difference may be the most direct factor in affecting gene expression, and it also indicated that there were no differences between the selected organs (late flower and early leaf bud) on gene expression responses to elevational changes. To some extent, all tissues and organs may cooperate in response to environmental changes.

The results of gene modules related to elevation based on a gene co-expression network analysis (WGCNA) showed negative correlation coefficients between modules and altitude (r < 0) ( Figure 4 ), indicating that genes in these modules were less expressed with higher altitude. We identified different environmental adaptation genes in the two organs. In the late flower bud, they were related to terpenoid metabolism ( Table 2 ; Figures 5A, C ), such as RsHMGR (DN1574_c1_g1) and RsTPS (DN13962_c0_g1, DN14863_c0_g1, etc.), which are key genes in terpenoid synthesis, while genes identified only in early leaf bud were mainly related to anthocyanin biosynthesis, leaf senescence and photosynthesis ( Table 2 ; Figures 5B, D ), such as the anthocyanin synthesis-related genes RsCHS (DN1880_c2_g1, DN31591_c0_g1, etc.) and RsF3’5’H (DN2888_c0_g1) ( Figure 6 ; Supplementary Figure S6 ). The expression levels of these genes decreased with increasing elevation. Terpenes are one of three main component classes of flower scent that is one of the most important traits for plant adaptation and evolution (Knudsen et al., 2006; Raguso, 2008; Wang et al., 2021a). Anthocyanins are secondary metabolites belonging to flavonoids determine the resistance of species to extreme environments (Grotewold, 2006; Tanaka et al., 2008). Environmental factors (temperature, light quality, humidity and oxygen concentration) change along an elevational gradient, with temperature and oxygen concentration decreasing with elevation, while light intensity (UV-radiation) increases (Körner, 2007). In Rhododendron plants, scent-related terpenes are the main components of their floral fragrance and the main substances that attract insects (Dudareva et al., 2005; Wang et al., 2021a). Wang et al. (2021a) demonstrated that species in Rhododendron sect. Azaleastrum, mainly distributed at low elevations, have fragrant flowers, but generally exhibit light or faded colors, compared to high-elevation rhododendrons, and the expansion of terpene synthase genes (TPSs) in R. ovatum not only promotes floral scent production but is also related to a greater diversity of the floral fragrant compounds. However, due to the limitation of abiotic factors, pollinators are generally rare in higher elevational environments (Hargreaves et al., 2015; Tong et al., 2021). Although the successful attraction of plant floral organs to pollinators determines the reproductive success rate of a species, to a certain extent, there must be an evolutionary trade-off for plants to adapt to the environment in this long evolutionary process (Armbruster, 2014). Because more fragrance production equates to more energy consumption, while it is extremely difficult for plant biomass accumulation where induced by harsh environmental factors at high elevations (Lavorel and Grigulis, 2012), when it is only to attract insects that are relatively rare there (Totland, 2001). It is also a balance for pay (energy expenditure on floral fragrance) and gain (attraction of pollinators). Notably, our study supported the hypothesis of a contrasting living strategy that were proposed for Trifolium repens (Hofmann and Jahufer, 2011) and Cyclocarya paliurus (Du et al., 2021). This is also a manifestation of plant adaptation to the environment as a long evolutionary process (Cruden, 1972; Hargreaves et al., 2015).

The accumulated anthocyanins are secondary metabolites of flavonoids (Tanaka et al., 2008; Jaakola, 2013), usually implicated in flower color as in the variety here (Ye et al., 2021). They are also frequently produced in young shoots, including leaf buds (Lee, 2002). Leaf buds often flush red synchronously during expansion that can be a widespread and visually striking phenomenon (Karageorgou and Manetas, 2006; Alberto et al., 2011). However, the red flushing appearance of leaf buds (as indicated in Figure 1A (3)) in many plants is transient, and most leaves contain considerable quantities of anthocyanins only in the juvenile stages (Dominy et al., 2002; Zhang et al., 2016), and can be a response to environmental stress (Lee, 2002). Among various environmental factors that affect anthocyanin biosynthesis, light intensity (UV radiation) and temperature are particularly important (Jaakola and Hohtola, 2010; Naing and Kim, 2021). It is generally acknowledged that strong light intensity can enhance the expression of genes related to anthocyanin biosynthesis and increase anthocyanin levels (Jaakola, 2013; Zhang et al., 2018). Studies have shown that within a certain range, light intensity can activate phytochromes and promote the synthesis or activation of a series of light-regulated enzymes such as phenylalanine lyase (PAL), chalcone synthase (CHS), and flavonoid glucosidase (UFGT), which affect the content and proportion of anthocyanins, thereby resulting in the coloration of leaves (Jaakola and Hohtola, 2010; Dong and Lin, 2021). However, the response of plants to light intensity varies in different environmental conditions. For example, Senica et al. (2017) revealed an increased anthocyanin content with increasing elevation in leaves of elderberry (Sambucus nigra), and Du et al. (2021) for Cyclocarya paliurus. Therefore, we hypothesize that increase on anthocyanin levels may be linked to the increase in light intensity across a certain elevational range. While at extremely high-elevation environments (> 3,000 m), gene expression is inhibited and hinders the synthesis and accumulation of anthocyanin in an early leaf bud. It has been shown that foliar anthocyanins are correlated with resistance to biotic and abiotic agents such as herbivores, cold and excess radiation (Gould et al., 2000; Close and Beadle, 2003). On the one hand, observations in the field showed that leaf buds, especially at lower elevation may be susceptible to herbivore damage, for insects may show preferences for young leaves for food or oviposition (Close and Beadle, 2003). Therefore, selective pressures of insect herbivory at lower elevations could elicit increased anthocyanin production (antiherbivore hypothesis) (Kost et al., 2020), because foliar anthocyanins and some other phenolics can serve as potent deterrents against generalist herbivores (Karageorgou and Manetas, 2006; Kost et al., 2020). On the other hand, excessive light intensities at high elevations may predispose leaves to photoinhibition of photosynthesis and may also cause damage to leaf buds because of their immature photosynthetic machinery (Miranda et al., 1981). However, Karageorgou and Manetas (2006) found that leaf thickness rather than anthocyanins could appreciably reduce the risk of photo-oxidative damage in young leaves in the field, and even the photoprotective function of anthocyanin accumulations in mature leaves under stress is transient. Thus, the effect of leaf thickness on photoprotection effectiveness may have been underestimated, which generally happened at high elevation. It could thus also be hypothesized that a substrate competition exists between lignin and anthocyanin biosynthesis (Hu et al., 2021). Plants occurring in high elevations can increase lignin deposits and thickness in the cell walls by competing for common substrates under UV exposure to reduce potential excessive light damage (Kost et al., 2020). In the present study, the down-regulation of DEUs involved in the anthocyanin metabolic pathway indicates that R. sanguineum var. haemaleum at high elevation experience a reduction in the expression of genes related to anthocyanin accumulation. This reduction is possibly induced by selective pressures of insect herbivory from environmental changes and substrate competition in the biosynthesis pathway. Our results are in accordance with a study on Zea mays where lowland maize landraces showed higher expression of several confirmed and putative genes involved in anthocyanin biosynthesis as compared to highland landraces (Kost et al., 2020). In addition, a high light intensity is also accompanied by low temperatures in high elevational environments, and studies have shown that low temperatures inhibit anthocyanin synthesis (Christie et al., 1994; Naing and Kim, 2021). Choi et al. (2009) reported that anthocyanin accumulation is light dependent at low temperatures and its synthesis is down-regulated at low temperatures. Collectively, the low expression level of anthocyanin biosynthesis-related genes at higher elevations can therefore be ascribed to a combined effect of biotic and abiotic factors.

The adaptation of organisms to elevational environments is a complex process (Abbott and Brennan, 2014; Sun et al., 2014). In plant species, some response mechanisms and related signalling systems are excited and reinforced to adapt to different environmental stresses (Weigel and Nordborg, 2015; Hao et al., 2019). As the two most important organs of plants, those of flowers (as reproductive organs) and leaves (as primary photosynthetic organs), representing the above-ground parts of the plant, their occurrence along an elevational gradient present, to a certain extent, the species’ responsiveness to environmental variation (Frei et al., 2014; Lopez-Goldar and Agrawal, 2021). In the present study, we found that the metabolic pathways and candidate genes related to altitude identified from flowers and leaves were significantly different ( Figure 5 ; Table 2 ). Genes from the late flower bud were mainly related to terpenoid metabolism, fatty acid elongation and flower organ development, while genes related to anthocyanin biosynthesis, photosynthetic antenna proteins and leaf senescence were identified in early leaf buds. Moreover, the expression of these genes was negatively correlated with altitude, that is, the higher the altitude, the lower their expression levels. When plants experience abiotic stress, they try to adapt to the stressful condition by adjusting the expression of a series of genes involved in complex networks (Nagano et al., 2019; Geng et al., 2021; Tokunaga et al., 2022).

Studies have shown that flower organ, as a very short-lived reproductive organ, mainly functions to attract pollinators, and produce seeds, while leaf organ is responsible for photosynthesis and energy supply (Bazzaz et al., 1979; Alderson and Rowland, 1995), and are long-lived in the evergreen Rhododendron sanguineum. Thus, it is not surprising that these two organ types behave very differently in response to changes in elevation. Terpenoids are the largest class of secondary metabolites in plants (Dudareva et al., 2005). Terpenes can be emitted from vegetative and floral organs to protect organs from heat burns produced by photosynthesis, and they can also participate in certain chemical ecological processes as signal molecules to attract insects and other pollinators (Pichersky and Gershenzon, 2002; Dudareva et al., 2004). Anthocyanins are secondary metabolites of flavonoids and play an important role in plant stress response (Tanaka et al., 2008; Du et al., 2021). The leaf is the main organ of plants exposed to air for photosynthesis, and it is more sensitive to the surrounding environment (including biotic and abiotic factors), mainly in the content of secondary metabolites (Hartmann, 1996; Landi et al., 2020). Moreover, under the long-term influence of external ecological factors, the morphological structure of leaves has great variability and plasticity, which is generally closely related to plant growth strategies and the ability to utilize resources (Obeso, 2002; Taylor et al., 2012; Rathee et al., 2021). Consequently, there are obvious functional divisions and constraints between the two organs from an ecological point of view.