Transcriptome analysis of Clematis lanuginosa: Novel features of the molecular events occurring under heat-shock stress

Clematis, known as the 'vine queen', is popular among flower lovers because of its large flowers, colorful and varied flower types and long flowering periods^[1,2]. The optimum growth temperature of the large-flowered group of Clematis is 15–25 °C. Heat damage occurs at greater than 30 °C, resulting in flower deformity, a shortened flowering period and withered branches^[3]. The study of heat-shock stress mechanisms in C. lanuginosa can increase our understanding of the heat-tolerance mechanisms in other large-flowered varieties of Clematis.

Transcriptome sequencing is a common method used in genetic research on various stresses encountered by plants, such as rice^[4], switchgrass^[5] and poplar^[6]. In Clematis apiifolia, a transcriptome analysis identified 1,074 up-regulated differentially expressed genes (DEGs) and 867 down-regulated DEGs after a heat-shock treatment. The heat-shock proteins, particularly small heat-shock proteins, are abundant under heat stress conditions^[7]. Under high temperature stress, 1,000–2,000 genes are up-regulated in heat-resistant rice germplasm^[8]. The regulatory mechanisms involved in the heat tolerance of plants are complex. When plants are exposed to high temperature stress, some genes are activated, resulting in increased levels of some metabolites and proteins. These genes may protect plants from high temperature stress, and some may activate the expression of certain genes to enhance the heat tolerance of plants. Heat-shock response is an important way to regulate heat tolerance through heat-shock factors (Hsfs) that regulate the accumulation of heat-shock proteins (Hsps)^[9]. The Hsps can be divided into Hsp100, Hsp90, Hsp70, Hsp60 and small heat-shock protein (sHsps) on the basis of their molecular weights^[10]. The sHSPs are chaperone proteins having molecular weights that range from 12 to 42 kD. Hsp gene expression is mainly regulated at the transcriptional level by Hsfs^[11].

Heat-stress granules, composed of cytosolic multi-chaperone complexes assembled by cytosolic sHSPs during heat stress^[12], become disintegrated during prolonged heat stress and recovery^[13]. When the formation of heat stress granules is blocked, plant resistance to high temperature stress is reduced^[13]. Ribosomal subunits are the main components of stress granules (SGs)^[14]. HSP90 facilitates the export of 60 S ribosomal subunits from resealed nuclear envelopes^[15], and HSP70 interacts with ribosomal subunits of thermotolerant lines in vitro^[16]. These results suggest that ribosomal subunits play important regulatory roles in heat tolerance.

In this research, we used transcriptome sequencing and a quantitative analysis to identify the genes of C. lanuginosa after 25 and 40 °C treatment. In total, six transcriptome databases were established and used to analyze differentially expressed genes under heat-shock stress conditions. The results provide a foundation for understanding the genetic mechanisms behind heat-shock stress responses in C. lanuginosa.

C. lanuginose seedings were used in this study. It was analyzed after 6 h at 25 and 40 °C. Under heat-shock conditions, the leaves slightly yellowed but did not dry up (Fig. 1a & b). The leaf relative water content decreased from 95.21% to 92.03% (Fig. 1c). The proline content increased from 32.98 µg/g to 38.57 µg/g (Fig. 1e). The NPQ value increased from 0.08 to 0.16 after heat-shock stress exposure (Fig. 1h). No significant differences were revealed in relative conductivity, Fv/Fm and Fv/Fo (Fig. 1d, f & g).

Figure 1. 

Photographs and physiological parameters of C. lanuginosa under heat-shock conditions. (a), (b) Photographs of C. lanuginosa under 25 °C and 40 °C. (c) Relative water content %. (d) Relative conductivity %. (e) Proline content(ug/g). (f) Fv/Fm. (g) Fv/Fo. (h) NPQ.

Library generation and clustering

To gain the C. lanuginosa. transcriptome under heat-shock stress conditions, three biological replicates each at 25 °C and 40 °C for 6 h were subjected to RNA-Seq. After pre-processing the sequence data of the six samples, 249,470,846 clean reads and 37,420,626,900 clean nucleotides (nt) were obtained in this experiment (Table 1). The Q20 percentages were between 96.35% and 97.15%. The GC percentages were between 46.38% and 47.96% (Table 1).

The high-quality reads were assembled into 91,090 contigs using the Trinity method. All the contigs produced 62,050 unigenes with an average length of 711 bp. The N50 length of the unigenes was 1,301 bp (Supplemental Table S1), and 16,439 unigenes were longer than 1 kb (Fig. 2). Six databases annotated 42,377 (68.27%) unigenes, NR (40,719, 65.60%), NT (30,237, 48.71%), Swiss-Prot (27,706, 44.64%), KEGG (25,033, 40.33%), COG (16,684, 26.88%) and GO (30,722, 49.50%) (Table 2).

Figure 2. 

Length distributions of contigs and unigenes in the transcriptome assembly.

GO enrichment and NR annotation

To functionally classify in all the unigenes, a Gene ontology (GO) enrichment of DEGs was analyzed using WEGO (Web Gene Ontology Annotation Plot) software. The GO annotated unigenes were assigned to cellular component (100,192) biological process (112,490) and molecular function (35,949) groups (Fig. 3), and 56 specific terms were significantly enriched (corrected P-value < 0.05). Most biological process genes were involved in 'Metabolic process' (20,106, 8.09%), 'Cellular process' (19,247, 7.74%) and 'Single-organism process' (12,530, 5.04%) (Fig. 3a). For the cellular component, 'Cell' (22,358, 8.99%), 'Cell part' (22,325, 8.98%) and 'Organelle' (10,807, 7.27%) were the dominant groups (Fig. 3b). For molecular function, 'Catalytic activity' (15,493, 6.23%), 'Binding' (14,799, 5.95%) and 'Transporter activity' (2,188, 0.88%) were the dominant groups (Fig. 3c).

Figure 3. 

GO classification analysis of unigenes in C. lanuginosa. The GO terms were summarized into three main categories: (a) biological process, (b) cellular component and (c) molecular function. The GO functions are shown on the left of the vertical axis, and the numbers and percentages of genes are shown on the right of the vertical axis.

The seven species with the most hits in the non-redundant (NR) annotation indicated that 78.47% of the unigenes' annotated sequences were from Solanum lycopersicum, Vitis vinifera, Prunus persica, Ricinus communis, Populus trichocarpa, Fragaria vesca subsp. Vesca, Glycine max, Cucumis sativus and Medicago truncatula. Among them, 38.19% of the unigenes were annotated from V. vinifera (Fig. 4).

Figure 4. 

NR protein sequence classification data. (a) E-value, (b) similarity, and (c) species distributions of the results of the NR annotation.

DEG enrichment analysis

To identify the DEGs associated with heat-shock stress in C. lanuginosa, six databases were constructed, three at 25 °C and three at 40 °C. In this research, 2,165 DEGs were found between the two treatments. Among these DEGs, 1,565 were up-regulated and 600 were down-regulated (Fig. 5).

Figure 5. 

Unigenes classified into expression classes. Red and blue genes are up-regulated and down-regulated, respectively.

The GO enrichment analysis revealed 'Cytoplasm' (778 DEGs), 'Cytoplasmic part' (676 DEGs), 'Biosynthetic process' (426 DEGs), 'Organic substance biosynthetic process' (406 DEGs), 'Cellular biosynthetic process' (398 DEGs), 'Cellular protein metabolic process' (342, DEGs), 'Organelle part' (346, DEGs) and 'Intracellular organelle part' (346, DEGs) (Fig. 6a). The results showed that the DEGs were primarily enriched in 'Cytoplasm', 'Cytoplasmic part', 'Biosynthetic process', 'Organic substance biosynthetic process', 'Cellular biosynthetic process', 'Cellular protein metabolic process', 'Organelle part' and 'Intracellular organelle part' pathways (Fig. 6a).

Figure 6. 

GO and KEGG analyses of DEGs in C. lanuginosa under heat-shock stress conditions. (a) GO analysis, (b) KEGG analysis.

The results of the KEGG enrichment analysis indicated that the DEGs were primarily enriched in 'Brite Hierarchies' (643 DEGs), 'Protein families: genetic information processing' (499 DEGs), 'Protein families: signaling and cellular processes' (183 DEGs) and 'Exosome' (118 DEGs) pathways (Fig. 6b).

qRT-PCR validation

To verify the DEGs under heat-shock stress conditions, we selected six important DEGs for a qRT-PCR analysis. These genes included Hsf30-like, Hsp70, sHsps and the Vacuolar protein sorting-associated protein gene (Supplemental Table S2). The qRT–PCR expression levels of these genes were in accordance with the transcriptomic data (Fig. 7). The results indicated that the transcription abundances of these genes under room temperature and heat-shock stress conditions were highly credible.

Figure 7. 

The expression levels of six genes in C. lanuginosa after 0 h and 6 h under heat-shock stress conditions using quantitative real-time PCR and RNA sequencing. Data are the means ± standard errors (n = 3).

Hsps and Hsfs in C. lanuginosa

In this research, Hsps and Hsfs were the most important gene families in heat-shock stress responses. In total, 51 unigenes were annotated as Hsps, including Hsp90, Hsp80, Hsp70 and sHsps. All the Hsps were up-regulated under heat-shock stress conditions. Among them, 35 were sHsps, accounting for 68.63% of the identified Hsps (Fig. 8a and Supplemental Table S3). Two Hsfs were identified as HsfA2 and Hsf24 (Fig. 8b and Supplemental Table S3), and they were significantly up-regulated under heat-shock stress conditions. One heat-shock cognate 70-kDa protein was also up-regulated in this research (Fig. 8c and Supplemental Table S3). These significantly regulated genes may play important roles in the heat-shock responses of C. lanuginosa.

Figure 8. 

Expression patterns of transcripts of heat stress-related genes. (a) Heat-shock proteins, (b) heat-shock factors, (c) heat-shock cognate 70-kDa protein 2, and (d) ribosomal proteins

Ribosomal proteins in C. lanuginosa

The expression levels of 12 ribosomal proteins were significantly up-regulated in this experiment (Fig. 8d and Supplemental Table S3). They were ubiquitin-40S ribosomal protein S27a, four 60S ribosomal protein L2 (Rpl2) and 30S ribosomal protein S (RpsG) 1, -2, -3, -4, -7, -11 and -13. The results indicated the ribosomal proteins may play important roles in the heat-shock stress regulation of C. lanuginosa.

Ethylene-responsive transcription factor expression pattern

Eleven ethylene-responsive transcription factors (Erfs) were found to be differentially expressed. Five, Erf4, -5, -053, Rap2-4 and Abr1, had up-regulated expression. Six genes (Erf1B, -010, -034, -053 and Rap2-12) had down-regulated expression (Fig. 9c and Supplemental Table S3).

Figure 9. 

Differentially expressed genes in C. lanuginosa. (a) Chalcone synthase, (b) Cysteine-rich receptor-like kinase, (c) Ethylene-responsive transcription factor and (d) Aspartic protease in guard cell.

Expression patterns of chalcone synthase, cysteine-rich receptor-like kinase and aspartic protease in guard cell

In our research, three chalcone synthase genes, chalcone synthase, -1 and -3, were detected as being significantly down-regulated (Fig. 9a and Supplemental Table S3). Eight cysteine-rich receptor-like kinases were differentially expressed under heat-shock conditions (Fig. 9b and Supplemental Table S3). Four unigenes annotated as aspartic protease in guard cell (Aspg1 and Aspg2) were down-regulated in our experiment (Fig. 9d and Supplemental Table S3).

The response mechanisms of plant heat tolerance are complex. On the basis of protein functions, the encoding genes can be divided into signal transduction (transcriptional genes and protein kinases) and functional (Hsps and catalases)^[17]. The corresponding signal mechanisms can be used to categorize the many heat tolerance-related signaling pathways in plants^[10]. The most characterized signal pathway is the Hsfs–Hsps pathway. The possible roles of Ca²⁺-dependency^[18], abscisic acid, salicylic acid^[19] and ethylene in heat-stress signaling have been studied. The Rps1-HsfA2-dependent chloroplast reverse regulatory mechanism is a new heat-stress regulatory mechanism^[20].

In our experiment, the Hsps gene family had the greatest up-regulated response. Hsps are chief functional proteins that are induced by heat-shock stress^[21]. Hsp expression is positively correlated with plant heat tolerance and improves the stress capacity of cells. In total, 34 Hsps have been identified as DEGs under heat-shock conditions in C. apiifolia, making it the most up-regulated gene family^[7]. Here, 51 Hsps were identified as DEGs in C. lanuginosa under heat-shock conditions, which made it the most up-regulated gene family. After the heat-shock treatment, the number of Hsps up-regulated in C. lanuginosa was greater than in C. apiifolia. Both C. apiifolia and C. lanuginosa are from heat resistant germplasm^[3]. However, in our follow-up observations, we determined that the heat resistance of C. apiifolia is stronger than that of C. lanuginosa. The Hsps represent a key gene family in heat- shock stress regulation, but the number of members may not be positively correlated with the plant heat tolerance.

SGs form when messenger ribonucleoproteins stall during translation initiation, and they can be induced by various cellular stresses, such as heat stress^[22,23]. The SGs contain a small ribosomal subunit (40S), translation initiation factors and specific RNA-binding proteins^[23]. In our research, 12 DEGs were annotated as ribosomal protein genes, including ubiquitin-40S ribosomal protein S27a, Rpl2 and Rps1, -4, -7, -11 and -13 (Fig. 8d). Mitochondrial HSP60s interact with WTF9 to regulate RNA splicing of ccmFC and rpl2 in Arabidopsis^[24]. In vitro chloroplast experiments have shown a more than two-fold increase of rpl23-rpl2 transcripts under heat-stress conditions^[25]. Rps13 is the most stable reference gene under high temperature conditions in Sesamia inferens^[26]. In this experiment, the differences among the expression levels of ribosomal protein genes after the heat-shock treatment indicated that these genes were involved in regulating the heat resistance of C. lanuginosa. The result indicated that Rps1-HsfA2-dependent chloroplasts reverse regulatory mechanisms exist and play important regulatory roles in heat-tolerance regulation in C. lanuginosa.

The Erf subfamily of transcription factors regulate the expression of genes related to salinity, drought, low temperature and biological stress, and they enhance stress tolerance in plants^[27−29]. Our research detected 11 ethylene-responsive transcription factors, indicating that this gene family participates in the regulation of heat-shock stress.

Chalcone synthase is a key enzyme in the flavonoid biosynthesis pathway^[30]. It plays important roles in the synthesis of plant secondary metabolites and in plant defense. Many stresses influence chalcone synthase expression, such as UV, pathogen attack and wounding^[30]. In recent years, studies have found that cysteine-rich receptor-like kinase is involved in plant disease resistance and defense responses^[31,32]. Aspartic proteinases in guard cell are a class of important proteolytic enzymes that play important regulatory roles in plant drought^[33]. These genes are rarely reported under heat-shock conditions. However, in this study, we found that they were affected by heat stress and differentially expressed, indicating that these gene-related expression pathways are involved in the regulation of heat tolerance of C. lanuginosa.

With the rise in global temperatures, it is urgent and necessary to breed heat-resistant plant varieties. Studying the heat-shock stress mechanisms of plants will lay the foundation for molecular breeding. The regulatory mechanisms of heat tolerance are complex. In this study, the parents of a large-flowered group of Clematis were selected as the research subject, and the regulatory mechanisms were analyzed. HsfA2 and Hsf24 were up-regulated Hsfs, sHsp were the up-regulated Hsps, and 12 ribosomal proteins may be involved in heat-shock stress regulation. The results lay a foundation for analyzing heat-resistant regulatory mechanisms and the molecular breeding of heat-resistant plants.

C. lanuginosa seedlings were collected from the Clematis germplasm nursery at the Institute of Botany, Jiangsu Province (China), and the Chinese Academy of Science, Nanjing, (China). Each seedling was pre-incubated in a light incubator at 25 °C for 7 d. Then, they were placed under either 40 °C (heat-shock stress) or 25 °C (control) conditions for 6 h. Three separate biological replicates were included for each treatment. The seedlings (approximately 10 g) were immediately frozen with liquid nitrogen, and then stored at –80 °C for subsequent analysis.subsequent analysis.

C. lanuginosa was identified by a known botanist Dr. Shuan Wang, who works at the Research Center of Ornamental Plants, Institute of Botany, Jiangsu Province and Chinese Academy of Science (China). The voucher specimen was deposited in the herbarium at the Institute of Botany, Jiangsu Province and Chinese Academy of Science (China). The certificate number is 000736328.break hyphen="true"/>cate number is 000736328.

RNA isolation and transcriptome library preparation

Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA). The total RNA content was determined using UV spectrophotometry, and the RNA quality was detected on 1% agarose gels. There were three biological replicates for each of the control and heat treatments. Finally, six separate paired-end databases were used for the RNA-sequencing (RNA-Seq) analysis. In accordance with the instructions, six paired-end libraries, with an average insertion length of 200 bp, were synthesized using a Genomic Sample Prep Kit (Illumina, San Diego, CA, USA). An Illumina HiSeq 2000 instrument was used to sequence libraries through a customer sequencing service (Novogene, Beijing, China).

RNA-Seq data processing, assembly and annotation

After RNA-Seq, raw reads were filtered and cleaned by eliminate low- quality raw reads and aptamer sequences. Non-redundant transcripts, comprising the clean reads, were assembled using Trinity software^[34]. Then, Q20, N percentage and GC percentage of the clean nucleotides were calculated. Clean, high-quality data were used for all subsequent analyses.

Parameter fragments per kilobase of transcript per million mapped fragments (FPKM) were used to quantitatively analyze gene expression levels.

Resulting sequences were used in BLAST searches and annotations against the National Center for Biotechnology Information database Clusters of Orthologous Groups of proteins (COG), nonredundant protein sequences (NR), Kyoto Encyclopedia of Gene and Genomes (KEGG) and protein sequence (Swiss-Prot) databases. Blast2go software was used for the gene ontology (GO) functional annotation. The raw data are available in the National Center for Biotechnology Information Sequence Read Archive (www.ncbi.nlm.nih.gov) under accession number PRJNA773622 (www.ncbi.nlm.nih.gov/bioproject/PRJNA773622).

Identification of differentially expressed genes

HTSeq v0.5.4p3^[35] was used to count read numbers mapped to each gene. On the basis of the gene length and the read counts, the reading per thousand bases was calculated. The P-values were adjusted using the Benjamin and Hochberg method to control the false discovery rate. The P values were adjusted using the Benjamini–Hochberg approach for controlling the FDR. Genes with an adjusted FDR ≤ 0.001 identified by DESeq^[36] and log2 FPKM (fold change) ≥ 1 were considered to be differentially expressed.

GO and KEGG enrichment analysis of differentially expressed genes

For the GO enrichment analysis, WEGO software^[37] was used to determine the statistical enrichment of the DEGs^[38,39]. The KEGG pathway analysis of DEGs was conducted using KOBAS software^[40] to determine the statistical enrichment of the DEGs. GO and KEGG analyses with corrected FDR ≤ 0.001 and log2 FPKM (fold change) ≥ 1 were considered significantly enriched DEGs.

qRT-PCR

We divided the data into two groups by comparing the data at 0 h with 6 h under 40 °C conditions. On the basis CDSs, we selected six unigenes with different expression levels for qRT-PCR verification. These unigenes were Unigene13812_Maoye, Unigene7274_Maoye, CL1187.Contig6_Maoye, CL4327.Contig1_Maoye, CL4528. Contig1_Maoye and Unigene20041_Maoye, and they were annotated as HSF30-like, HSP68, 18.2-kDa, 16.9-kDa and 15.7-kDa class I sHsps, and Vacuolar protein sorting-associated protein 28, respectively (Supplemental Table S2). The primer sequences are shown in Supplemental Table S4. This primers were designed using Beacon Designer 8.0 software in accordance with the qRT-PCR primer requirements. A StepOne instrument (Applied Biosystems Company, USA) was used in this experiment. The amplification procedure was as follows: step 1: 95 °C for 30 s; step 2: 40 cycles of 95 °C for 5 s and annealing at 60 °C for 30 s; step 3: 95 °C for 15 s, 60 °C for 60 s and then 95 °C for 15 s. Each qRT-PCR was performed with three biological repeats.