Transcriptome analysis reveals key genes in response to high-temperature stress in Rhododendron molle

Abstract

Rhododendron molle is a deciduous rhododendron, a high-altitude plant prized for its medicinal and ornamental properties. A major challenge when introducing this plant to lower altitudes is understanding its response to high-temperature stress. Using transcriptome analysis, this study examined leaves under varying temperatures, identifying 344,593 transcripts, 124,901 Unigenes, and 12,089 differentially expressed genes (DEGs) at 36 °C high-temperature stress (ST). At 42 °C high-temperature stress (SY), 12,032 DEGs were found, indicating a significant impact of temperature on gene expression. A Gene Ontology analysis (GO) revealed that these DEGs are mostly involved in stress response, catalytic activity, binding, transporter activity, and immune processes. A Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis highlighted enrichment in pathways such as plant-pathogen interaction and spliceosomes, suggesting their key roles in the temperature stress response. Key genes such as Brassinosteroid-Insensitive 1-Associated Receptor Kinase (BAK1) and Pathogenesis-Related Gene Transcription Activation Factor (PTI6) were upregulated under ST stress, whereas heat shock proteins (HSP83A) and HSP90-1 were downregulated under SY. These findings offer insights into the molecular response of Rhododendron molle to high temperature, aiding further research in this area and potentially improving the plant’s cultivation and application.

Full Article

Transcriptome Analysis Reveals Key Genes in Response to High-Temperature Stress in Rhododendron molle

Linshi Wu,^a,b,* Yan Liu,^a,b,* Xinyun Liu,^a,b Qiaoyun Li,^a,b Xinyu Yi,^a,b Chan Chen,^a,b Ling Wang,^a,b and Juyang Liao ^a,b

Rhododendron molle is a deciduous rhododendron, a high-altitude plant prized for its medicinal and ornamental properties. A major challenge when introducing this plant to lower altitudes is understanding its response to high-temperature stress. Using transcriptome analysis, this study examined leaves under varying temperatures, identifying 344,593 transcripts, 124,901 Unigenes, and 12,089 differentially expressed genes (DEGs) at 36 °C high-temperature stress (ST). At 42 °C high-temperature stress (SY), 12,032 DEGs were found, indicating a significant impact of temperature on gene expression. A Gene Ontology analysis (GO) revealed that these DEGs are mostly involved in stress response, catalytic activity, binding, transporter activity, and immune processes. A Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis highlighted enrichment in pathways such as plant-pathogen interaction and spliceosomes, suggesting their key roles in the temperature stress response. Key genes such as Brassinosteroid-Insensitive 1-Associated Receptor Kinase (BAK1) and Pathogenesis-Related Gene Transcription Activation Factor (PTI6) were upregulated under ST stress, whereas heat shock proteins (HSP83A) and HSP90-1 were downregulated under SY. These findings offer insights into the molecular response of Rhododendron molle to high temperature, aiding further research in this area and potentially improving the plant’s cultivation and application.

DOI: 10.15376/biores.19.4.8238-8256

Keywords: Rhododendron molle; High temperature stress; Transcriptome analysis; Differentially expressed genes (DEGs)

Contact information: a: Hunan Botanical Garden, Box 4100004, Changsha P.R. China; b: Rhododendron Engineering Technology Research Center, Box 4100004, Changsha P.R. China; *These authors contributed equally to this work; Corresponding authors: Juyang Liao liaojuyang@163.com (Hunan Botanical Garden, Changsha, PR China)

INTRODUCTION

Rhododendron molle, distinctively known as the “yellow azalea” or “noisy sheep flower,” holds the singular distinction of being the only species native to China within the subgenus Pentanthera of the Rhododendron family (Ericaceae). Celebrated for its exceptional ornamental appeal, it plays a crucial role in enriching the color diversity of azaleas (Ureshino et al. 1998). This plant predominantly flourishes at altitudes of around one thousand meters, with its main habitat spread across central and southern China and a smaller presence in the southwestern region. Notably, Rhododendron molle is currently recognized as an endangered species (Wang et al. 2010), underscoring the need for conservation efforts to protect this unique and valuable botanical treasure.

Rhododendron molle, a plant with significant medicinal qualities, is enriched with compounds such as lignans in its flowers and flavonoids, rhododendron toxins, and methyl-lichenic acid in its leaves. Lignans are an important kind of natural phytoestrogen with various biological activities, including antioxidant, antiviral, and anti-tumor activities. The lignans in rhododendrons have little effect on animals, but the gray anotoxins have a strong insecticidal effect and have been widely used as a pesticide. When applied via contact or fumigation, they poison the stomach, interfere with feeding, and inhibit the growth and development of many kinds of pests. Rhododendron molle’s other constituents endow it with diverse therapeutic benefits. Traditionally, it has been used for its effectiveness in dispelling wind, alleviating pain, and reducing dampness, and it also serves as an anesthetic. It is particularly effective in treating conditions such as fractures, rheumatism, and persistent skin eczema (Liu et al. 2007; Cheng et al. 2011; Zhi et al. 2013; Li et al. 2018). Intriguingly, certain toxins in Rhododendron molle, especially those derived from the “noisy sheep flower,” have shown remarkable agricultural value. They demonstrate potent stomach-poisoning, contact-killing, and fumigation effects against various pests (Zhu et al. 2023). Notably, research has highlighted that these toxins are effective in deterring feeding in the larvae of pests like Leptinotarsa decemlineata and Spodoptera frugiperda, underscoring the plant’s potential in natural pest control (Shin and Yu 1993). Aside from its medicinal properties, Rhododendron molle also stands as a pivotal resource in azalea breeding. Azalea breeders have long recognized its excellent compatibility with other azalea species, making it an ideal candidate for hybridization (Geng et al. 2019). Many of its traits, typically recessive in the first generation, render it a superior parent for breeding programs. Despite its critical role in azalea cultivation, the challenges presented by Rhododendron molle such as its low seed germination rate and the difficult rooting and low survival rates of cuttings, have hindered its use in urban landscaping (Meijón et al. 2009). Consequently, intensifying research to improve Rhododendron molle cultivation techniques and promote their broader implementation is crucial for diversifying urban tree species and holds significant implications for agricultural production and medical research (Zhou and Zhu 2020).

Rhododendron molle grows under shrubs or mixed forests on mountain ridges in mountainous grasslands or hilly areas at an altitude of 1000 meters, with the most suitable growth temperature being 15 to 20 ℃; these conditions adversely affect the introduction and garden cultivation of Rhododendron molle (Arisumi et al. 1989). Additionally, global climate change has led to more frequent occurrences of hot weather and water shortages across the nation, impacting the normal growth of plants. This situation of high temperatures is expected to worsen in the future (Schlenker and Maximilian 2018; Yuan et al. 2019). Therefore, it is crucial to understand and explore the molecular regulatory mechanisms that enable Rhododendron molle to withstand high-temperature stress, which is particularly important for breeding purposes.

RNA-Seq can be used to comprehensively obtain the sequence and expression information of the transcriptome of a particular organism under different environmental conditions and analyze the differential expression of genes at the RNA level to identify key DEGs. This, in turn, allows an in-depth study of the regulatory mechanisms and functions of differential genes in different environments (Yan et al. 2004). For example, 18 genes encoding heat stress proteins (HSPs) were identified in Paeonia suffruticosa by means of transcriptome analysis, and all of them showed significant up-regulation under high-temperature treatment. A transcriptome analysis of Narcissus tazetta var. chinensis in response to high-temperature stress identified 96 known miRNAs, 25 novel miRNAs and 52,286 predicted siRNAs that may play important roles in response to heat stress in daffodils (Yan et al. 2004). In this study, we analyzed the transcripts of Rhododendron molle under high-temperature stress using RNA-seq technology to obtain the differentially expressed genes. Our results can provide a reference for the later molecular breeding of Rhododendron molle variants that are robust against high temperatures.

EXPERIMENTAL

Study Sites and Plant Materials

The experimental materials were kept in Hunan Provincial Botanical Garden (113°01′30″E, 28°06′40″N) Two-year-old cuttings of Rhododendron truncatulum were transplanted into pots (h≥20cm), with one cutting planted in each pot. A total of 2 high-temperature treatments (36 ℃ and 42 ℃ referred to as ST and SY, respectively) were set up, and a 25 ℃ treatment was used as the control (CK). The experimental design included lighting at 2000 lx(14 h day/10 h night) in an artificial climate chamber. Three biological replicates were planted in the control and each treatment. Under continuous high-temperature stress in the SY treatment, Rhododendron molle displayed symptoms of yellowing and blackening of its leaves after 10 days, with signs of death observed after 14 days. Therefore, a high-temperature stress duration of 15 days was designed (Jiang et al. 2023). On the 15th day, healthy leaves without pests or diseases were collected from the target plants, rapidly frozen in liquid nitrogen, and stored at -80 °C. These samples were then sent to Nanjing Jisi Huiyuan Biotechnology Co., Ltd. for further analysis.

RNA Extraction and Sequencing Library Construction

The mRNA was extracted from the biological samples using the Poly (A) structure, and the RNA quality was determined using a Nanodrop ultra-micro UV-Vis spectrophotometer and Agilent 2100 Bioanalyzer. Eukaryotic mRNA was enriched using magnetic beads with Oligo (dT), and the first cDNA strand was synthesized by adding Fragmentation Buffer to randomly interrupt the mRNA, followed by the synthesis of the double-stranded cDNA strand, purification of the cDNA, and end-repair, etc. Finally, the cDNA library was enriched via PCR.

Assembly and Functional Annotation of the Transcriptome

High-throughput sequencing (e.g., Illumina HiSeq 2500) was used for experimental processing. The Unigene sequence were compared with the NR (Marchler et al. 2015), Swiss-Prot (Gasteiger et al. 2001), GO (Gene Ontology Consortium, 2004), COG (Tatusov et al. 2003), KOG (Jensen et al. 2007) and KEGG (Kanehisa et al. 2012) databases using BLAST v2.3.0+ (http://blast.ncbi.nlm.nih.gov/Blast.cgi) software, and the Unigene amino acid sequence were predicted and compared with the Pfam database (Mistry et al. 2021) using HMMER v3.4 (Finn et al. 2011) software to obtain the Unigene annotation information.

Functional Analysis of the DEGs

The sequenced reads were compared with the Unigene library using Bowtie2 v2.4.4 (Langmead and Salzberg 2012), and the expression level was estimated based on the comparison results in combination with RSEM (Li and Dewey 2011). The expressing abundance of the corresponding Unigene is expressed using FPKM values.

RESULTS AND DISCUSSION

Sequencing Data Assembly Results

Each sequencing sample was repeated once, and after sequencing quality control, the accuracy of each sample base was greater than 96% (corresponding to Q20 ), and the GC content was basically below 50%. Thus the sequencing data could be used for assembly (Moreton et al. 2016). The data output and assembly statistics for each sample are shown in Table 1. For the de novo assembly of the transcriptome, Trinity assembly software was adopted. The assembly resulted in a total of 124,901 Unigenes, with a total length of 530,057,171 nt. and an average length of 1538.21 nt, N50. All the Unigenes were arranged in order from longest to shortest, and their lengths were added. When a fragment’s length is 50% of the total number of fragments, the length and number of the fragment, i.e., the length and number of N50, indicate the quality of the assembly: the lower these values are for the Unigenes, the better the quality.

Table 1. Sequencing Sata Quality

Transcriptome Unigene Cluster Results

Since there were no reference genomic data for Rhododendron molle, Clean reads were used in Trinity software for hybrid splicing. A total of 344,593 transcripts and 124,901 Unigenes were obtained. The average length of the transcripts was 1538.21, and the N50 values of the transcripts and Unigenes were 2466 and 1725, respectively (Table 2), indicating high assembly integrity. This data met the requirements for the next step of gene function annotation and classification.

Table 2. Data Assembly Statistics

The assembled transcriptome sequences were compared with COG, GO, KEGG, KOG, Swiss-Prot, Pfam, NR, and other databases to obtain the annotation information of each database. The numbers of Unigenes annotated in the above databases were 10,973, 32,259, 51,672, 29,049, 23,701, 35,661, 62,917, respectively. The gene function annotation statistics are shown in Table 3.

Table 3. Annotation Statistics of Single Gene Clusters

Analysis of Differentially Expressed Genes Under Different Heat Stresses

Based on the RPKM values, the transcript gene expression differences under the different treatments were assessed, with screening criteria of FDR≤0.05 and FC≥2. In Fig. 1, the horizontal coordinate indicates the logarithmic value of the difference multiplier between the two treatment groups, and the vertical coordinate indicates the negative log10 value of the FDR of the difference between the two subgroups. In these volcano plots of differentially expressed genes, red (up-regulated expression) and green (down-regulated expression) dots indicate differential gene expression, while gray dots indicate no difference. By comparing the differential expression of genes between high-temperature stress and the control, the volcano plot analysis (Fig. 1) showed a total of 12,089 differentially expressed genes, of which 6004 were up-regulated and 6,085 were down-regulated, in the control compared with the 36 °C high-temperature stress. Additionally there were 12,032 differentially expressed genes in the control compared to the 42 °C high-temperature stress, of which 5,431 were up-regulated and 6,601 were down-regulated. There were 4594 common differentially expressed genes between the two groups of differential genes (Fig. 2), of which 1224 were co-up-regulated and 3007 were co-down-regulated, in addition to 363 common reverse differentially expressed genes between the two. Some of the common differentially expressed genes responded to both 36 °C and 42 °C high-temperature stress. The number of down-regulated genes in Rhododendron molle was higher than the number of up-regulated genes with increasing temperature under the high-temperature stress treatments, indicating that high-temperature stress significantly suppressed some of the physiological activities of Rhododendron molle.

Fig. 1. Volcano plot of the number of differentially expressed genes. Each point in the differential expression volcano plot represents a gene, and the horizontal coordinate indicates the logarithmic value of the expression ploidy of a gene in the two samples; the vertical coordinate indicates the negative logarithmic value of FDR. The larger the absolute value of the horizontal coordinate, the greater the difference in expression ploidy between the two samples; the larger the value of the vertical coordinate, the more significant the differential expression and the more reliable the differentially expressed genes obtained from the screening. The green dots represent down-regulated differentially expressed genes, the red dots represent up-regulated differentially expressed genes, and the gray dots represent non-differentially expressed genes.

Fig. 2. Venn diagram of the number of differentially expressed genes

GO Classification and Enrichment Analysis of Significantly Different Genes

The obtained target genes were subjected to GO enrichment analysis. The biological functions of the differentially expressed genes were mainly analyzed in terms of their molecular functions, cellular components, and biological processes. The biological processes enriched on metabolic processes, cellular processes, stimulus responses, and biological regulation; the cellular component part enriched on cells, cellular components, organelles, organelle components, membranes, membrane components, and protein complexes; and the molecular function part enriched on binders and catalytic activity. The GO enrichment results for Rhododendron molle under different high-temperature stresses are displayed in Fig. 3. The significantly enriched GO term categories did not change, but the numbers of enriched genes all increased, indicating that high-temperature stress led to increased plant cell destruction.

Fig. 3. Differential gene GO functions. Cellular component 1: Extracellular region; 2: Cell; 3: Nucleoid; 4: Membrane; 5: Cell junction; 6: Membrane-enclosed lumen; 7: Protein-containing complex; 8: Organelle; 9: Other organism; 10: Other organism part; 11: Extracellular region part; 12: Organelle part; 13: Membrane part; 14: Cell part; 15: Supramolecular complex . molecular function 1: Catalytic activity; 2: Structural molecule activity; 3: Transporter activity; 4: Binding; 5: Antioxidants activity; 6: Nutrient reservoir activity; 7: molecular transducer activity; 8: molecular function activity; 9: molecular carrier activity; 10: Transcription regulator activity. Biological process: 1: Reproduction; 2: Cell killing; 3: Immune system process; 4: Metabolic process; 5: Cell proliferation; 6: Cellular process; 7: Carbon utilization; 8: Nitrogen utilization; 9: Reproductive process; 10: Biological adhesion; 11: Signaling; 12: Multicellular organismal process; 13: Developmental process；14: Growth；15: Locomotion；16: Pigmentation；17: Rhythmic process; 18: Response to stimulus; 19: Localization; 20: multi-organism process; 21: Biological regulation; 22: Cellular component organization; 23: Detoxification.

KEGG Enrichment Analysis of Significantly Different Genes

The KEGG database was used to analyze the enrichment pathways of the differentially expressed genes detected in Rhododendron molle under different high-temperature stress treatments. As a result 4594 differentially expressed genes (DEGs) were found to be enriched in 125 pathways. The KEGG metabolic pathways for the control in comparison to 36 °C high-temperature stress mainly involved sulfur metabolism, astragalus, diarylheptane and gingerol biosynthesis, starch and sucrose metabolism, plant-pathogen interactions, phytohormone signaling, photosynthesis tentacle protein, phenylpropanoid biosynthesis, and other pathways. Among the top 20 metabolic pathways with the most reliable enrichment significance, as shown in Table 4, the plant-pathogen interaction pathway was annotated to a maximum of 88 genes; followed by spliceosomes (annotated to 60 differentially expressed genes), plant hormone signaling (annotated to 59 differentially expressed genes), ribosomes (annotated to 58 differentially expressed genes), and other pathways. The KEGG metabolic pathways for the control in comparison to 42 °C high-temperature stress mainly involved in plant-pathogen interactions, ribosomes, spliceosomes, endoplasmic reticulum protein processing, and other pathways. Among the top 20 metabolic pathways with the most reliable enrichment significance, the plant-pathogen interaction pathway was annotated to the most genes (100), followed by the ribosomal (annotated to 83 differentially expressed genes), ribosomal (annotated to 73 differentially expressed genes), and endoplasmic reticulum protein processing (annotated to 73 differentially expressed genes) pathways.

Table 4. KEGG Pathway Enrichment Information for the Top 10 Differentially Expressed Genes

Notes: ST: 36 °C high-temperature stress; SY: 42 °C high-temperature stress

Fig. 4. KEGG functional enrichment analysis of differentially genes

Functional Identification of Differentially Expressed Genes

From the enrichment results for differentially expressed genes in the KEGG signaling pathways, the related genes that were up-regulated or down-regulated in all metabolic pathways of lambda R. molle under high-temperature stress were identified. Two significantly up-regulated genes were screened from among the most enriched genes in the plant-pathogen interactions pathway under the ST high-temperature treatment(36 °C): oleuropein receptor protein-associated receptor kinase (BAK1) and disease process-associated gene transcriptional activator (PTI6).Four genes were significantly down-regulated, including mitogen-activated protein kinase kinase 1a (MKK1a), mitogen-activated protein kinase kinase kinase 5 (MKK5), and retroviral-associated Pol polyprotein (RE1) (Tables 5 through 7). Screening under the spliceosome pathway showed significant upregulation of the sulfate glycoprotein (RH14) and RNA unwinding helicase (SMD1A) and significant downregulation of a small nuclear ribonucleoprotein (STA1). Under the phytohormone signal transduction pathway, indole-3-acetic acid amide synthase (GH3.10), two-component response regulator (RR22), and DELLA protein (GAI1) were the three up-regulated genes and gene erythromycin receptor (GID1B) and oleuropein receptor protein-associated receptor kinase (BAK1) were the two down-regulated genes. Oleuropein receptor protein-related receptor kinase (BAK1), WRKY transcription factor (WRKY22), and transcriptional activator of disease process-related genes (PTI6) were significantly up-regulated in the phytopathogen interactions pathway under the SY high-temperature treatment (42 °C), while heat shock protein HSP83A, heat kinin (HSP90-1), and mitogen-activated protein kinase kinase 5 (MKK5) were significantly downregulated. The ribosomal proteins RPS3C and RPL24 showed up-regulation while RPS20B RPL10A and showed down-regulation under the ribosomal pathway. Protein homologs (MAGO1) showed up-regulation and serine-arginine-rich splicing factors (SC35, SR34, SR34B, SCL33) all showed down-regulation under the spliceosomal pathway (Tables 8-10). Genes in the metabolic pathway exhibited different levels of coordination involving multiple pathways such as photosynthesis, respiration, redox, phytohormone signaling, and gluconeogenesis, which work together to combat high-temperature stress.

Table 5. Information on Relevant DEGs Under the Interactions Pathway in Plants

Table 6. Information on Relevant DEGs Under the Spliceosome Pathway

Table 7. Information on Relevant DEGs Under the Phytohormone Signaling Pathway

Table 8. Information on Relevant DEGs Under the Plant-Pathogens Interactions Pathway