Jiahong Xu1,2,3
Meng Fang1,2,4Zhihao Li1,2,3Maoning Zhang1,2,3Xiaoyu Liu1,2,3Yuanyuan Peng1,2,3
Yinglang Wan3
Jinhui Chen1,2,3*- 1Key Laboratory of Genetics and Germplasm Innovation of Tropical Special Forest Trees and Ornamental Plants, Ministry of Education/Engineering Research Center of Rare and Precious Tree Species in Hainan Province, College of Forestry, Hainan University, Haikou, China
- 2Hainan Key Laboratory for Biology of Tropical Ornamental Plant Germplasm, College of Forestry, Institute of Tropical Agriculture and Forestry, Hainan University, Haikou, China
- 3Hainan Key Laboratory for Sustainable Utilization of Tropical Bioresources, College of Tropical Crops, Hainan University, Haikou, China
- 4School of Life Sciences, Lanzhou University, Lanzhou, China
Long non-coding RNAs (lncRNAs) regulate plant responses to abiotic stresses. However, the short reads produced by second-generation sequencing technology make it difficult to accurately explore full-length transcripts, limiting the study of lncRNAs. In this study, we used third-generation long-read sequencing technology with the PacBio Sequel and Illumina platform to explore the role of lncRNAs in the heat stress response of Populus x canadensis Moench trees. We using 382,034,416 short reads to correct 4,297,179 long reads by resulted in 66,657 full-length transcripts, representing 33,840 genes. Then, 753 putative lncRNAs were identified, including 658 sense lncRNAs (87.38%), 41 long intervening/intergenic non-coding RNAs (lincRNAs) (5.44%), 12 antisense lncRNAs (1.59%), and 42 sense intronic lncRNAs (5.58%). Using the criteria | log2FC| ≥ 1 and q-value < 0.05, 3,493 genes and 78 lncRNAs were differentially expressed under the heat treatment. Furthermore, 923 genes were detected as targets of 43 differently expressed lncRNAs by cis regulation. Functional annotation demonstrated that these target genes were related to unfolded protein binding, response to stress, protein folding, and response to stimulus. Lastly, we identified a lncRNA–gene interaction network consisting of four lncRNAs and six genes [Heat Shock Protein 82 (HSP82), HSP83, Disease Resistance Protein 27 (DRL27), DnaJ family protein (DNJH), and two other predicted protein-coding genes], which showed that lncRNAs could regulate HSP family genes in response to heat stress in Populus. Therefore, our third-generation sequencing has improved the description of the P. canadensis transcriptome. The potential lncRNAs and HSP family genes identified here present a genetic resource to improve our understanding of the heat-adaptation mechanisms of trees.
Introduction
High temperature is a common abiotic stress experienced by plants. Approximately 23% of land on the Earth has an annual mean air temperature above 40°C (Leone et al., 2003). The frequent high temperatures that have occurred throughout the world in recent years have seriously threatened the growth of plants (Mora et al., 2015). High temperatures alter plant metabolic processes, thus affecting growth and development. Plant adaptations to acclimatize to high temperature include changes in morphology and physiology. For example, in response to heat stress, plants have shown decreased photosynthesis and transpiration (Mathur et al., 2014), and cell membrane stability (Havaux et al., 1996) and increased levels of antioxidants (Xu et al., 2006; Zou et al., 2011), phytohormones (Abdelrahman et al., 2017), compounds used for osmotic adjustment (Wahid et al., 2007), and heat shock proteins (HSPs) (Miller and Mittler, 2006). The changes at the molecular level take place immediately after the onset of heat stress and result in altered gene expression and transcript accumulation, and the synthesis of HSPs to induce stress tolerance (Schulze et al., 2005). HSPs enhance protein stability by helping other proteins to not unfold during heat stress (Wahid and Close, 2007).
HSPs are mainly divided into six classes based on their molecular weights: HSP100, HSP90, HSP70, HSP60, HSP40, and small HSPs (Kotak et al., 2007). In plants, the expression of many HSP genes increases and regulates the synthesis of HSPs after exposure to heat stress (Ada et al., 1987; Kant et al., 2008). For example, previous transcriptome analysis studies reported similar changes in expression of HSP genes in Lolium perenne and Festuca arundinacea (Hu et al., 2014; Wang K. et al., 2016). HSP40 and HSP70 were up-regulated in Dianthus caryophyllus under heat stress (Wan et al., 2015). In addition, transgenic creeping bentgrass (Agrostis stolonifera) plants expressing a rice (Oryza sativa) SUMO ligase showed increased expression of two HSP genes and improved resistance to high temperatures compared with non-transgenic plants (Li et al., 2013). These results suggested that HSPs play an important role in response to high temperature stress in plants.
LncRNAs have many key regulatory functions in plants. For example, a cold-inducible intronic lncRNA silences the floral repressor FLOWERING LOCUS C (Heo and Sung, 2011). Moreover, lncRNAs have roles in abiotic stress responses. For example, overexpression of the lncRNA npc536 in Arabidopsis thaliana enhances plant salt tolerance (Amor et al., 2009). In addition, the targets of rice stress-induced lncRNAs were enriched in genes involved in development and stress tolerance, suggesting that lncRNAs also play important roles in rice growth under stress (Yuan et al., 2018).
Upon heat stress, lncRNAs play an important role. For example, lncRNA HSR1 associates with protein eEF1A in response to heat stress (Shamovsky and Nudler, 2009). Nucleolar sequestration of HSP72 is facilitated by lncRNA in response to heat stress (Kotoglou et al., 2009). In addition, Transcription elongation of HSPs were regulated by lncRNA, and suppression of transcription after heat stress were mediated by lncRNA (Yakovchuk et al., 2009, 2011; Place and Noonan, 2014). Although there is a lot of research on lncRNA, our understanding of the role of lncRNA in heat stress is very limited. Development of research methods and techniques, especially for the identification and characterization of lncRNAs, will improve our understanding of non-coding RNA functions. Many studies have reported on the response of regulatory lncRNA genes to stress using RNA sequencing (RNA-seq), which enables accurate, genome-wide analysis of transcriptomes (Au et al., 2012). However, high-throughput sequencing produces short reads; the assembly of transcripts from these reads can be substantially improved using the long reads produced by single-molecule long-read “third-generation” sequencing technologies. Indeed, many studies have harnessed third-generation sequencing (Treutlein et al., 2014; Abdel-Ghany et al., 2016; Bo et al., 2016) and found that the resulting full-length sequences are more accurate for exploring gene length and lncRNAs (Roberts et al., 2013; Hackl et al., 2014).
Here, we used single-molecule real-time (SMRT) technology combined with high-throughput sequencing to reveal changes in lncRNA abundance in Populus x canadensis Moench after treatment at 40°C for 1 h. To the best of our knowledge, this study using third-generation sequencing to characterize the P. canadensis transcriptome and serves as a basis for further research on the response to high temperature.
Materials and Methods
Plant Materials and Heat Stress Treatment
For the heat stress treatment, we selected uniformly developed seedlings (one month) of P. canadensis and transplanted them into individual pots. The seedlings were grown in a greenhouse at Hainan University (Danzhou; 109°29′ 25″ E, 19°30’ 40″ N) and grown for five months, and six uniform individuals were used for the treatment. The plants were named pch1, pch2, and pch3 (for P. canadensis heat) for the heat-stress treatment under 40°C for 1 h, and pcq1, pcq2, and pcq3 (for P. canadensis control) for the untreated controls. Leaf samples were collected from the six plants for RNA extraction. The collected leaves were immediately frozen in liquid nitrogen and stored in a −80°C ultra-low temperature freezer until they were used for RNA isolation.
RNA Isolation and Illumina Sequencing
Total RNA of mature leaves collected from pcq and pch groups was extracted by CTAB method. The quality of the six RNA samples was assessed by NanoDrop 2000 (Thermo Scientific) and Agilent 2100 Bioanalyzer (Agilent Technologies). Only RNA samples with absorption values in the 260/280 ratio from 1.9 to 2.2, and in the 260/230 ratio ≥2.0, and we used only samples that had an RNA integrity number (RIN) >6.8 for follow-up experiments. Polyadenylated mRNA was enriched by oligo (dT) magnetic beads.
For the short-read RNA-seq, fragmentation buffer was added to fragment the mRNA into shorter pieces. We synthesized single-stranded cDNA from the mRNA using random hexamer primers and synthesized double-stranded cDNAs by adding buffer, dNTPs, and DNA polymerase I. We purified the double-stranded cDNA using AMPure XP beads and subjected it to end repair, addition of the poly-A tail, ligation of the sequencing linker, and fragment size selection. Finally, the cDNA libraries were subjected to PCR enrichment and sequenced with the Illumina HiSeq 2500.
PacBio Iso-Seq Library Preparation
For SMRTbell libraries, were combined equal amounts of total RNA from the biological replicates and generated two RNA pools, the heat stress pool and the control pool. From these pools, oligo (dT) was used to enrich for mRNA containing a poly-A tail, then mRNA was reverse transcribed into cDNA using the SMARTer PCR cDNA Synthesis Kit. We used PCR to amplify the cDNAs. The fragments were then screened for large-scale PCR to obtain sufficient cDNA. The full-length cDNA obtained was subjected to injury repair, end-repair, ligation of SMRT dumbbell-type linkers, and construction of a full-length transcriptome library. We removed the sequence of the unligated linker at both ends of the cDNA and then bound the primer and the binding DNA polymerase to form a complete SMRTbell library.
The libraries were then sequenced using PacBio Sequel II System and SMRT. The raw Iso-Seq data were processed on SMRTlink v6.0 software to obtain subread sequences. A circular consistency sequence (CCS) was obtained from the correction between subreads. The full-length sequences that contained a 5′ primer, 3′ primer, and poly-A tail were clustered using the Iterative Isomer Clustering (ICE) algorithm. Finally, the obtained consensus sequence was calibrated using the short-read sequences to obtain a high-quality sequence for subsequent analysis.
Identification and Functional Annotation of LncRNAs
To investigate the function of genes, we used BLAST (version 2.2.26) (Altschul et al., 1997), HMMER 3.1 (Eddy, 1998), and KOBAS 3.0 (Xie et al., 2011) to search the NR (Li et al., 2002), Nt (Li et al., 2002), KOG (Koonin et al., 2004), COG (Tatusov et al., 2003), Swiss-prot (Amos and Rolf, 2000), KEGG (Minoru et al., 2004), GO (Ashburner et al., 2000), and protein family (Pfam) (Finn et al., 2016) databases. CPC v0.9 (Kong et al., 2007), CNCI v2 (Sun et al., 2013), PLEK v1.2 (Aimin et al. (2014), and Pfam-scan (Finn et al., 2016) were used to identify lncRNAs from the Iso-Seq data (Chao et al., 2018). We use CNCI with default parameters. CPC (set the e-value “1e-10”) was used to assess the extent and quality of the ORF in a transcript and search the sequences with NCBI eukaryotes’ protein database to clarify the coding and non-coding transcripts. Pfam searches use default parameters of −E 0.001 −domE 0.001. The PLEK used to assess coding potential for species that lack high-quality genomic sequences and annotations. PLEK used default parameters of -minlength 200, and removed transcripts <200 bp. Transcripts predicted with coding potential by either/all of the three tools above were filtered out, and those without coding potential were our candidate set of lncRNAs. Since the whole genome sequencing of P. canadensis has not yet been completed, lncRNAs were classified into four categories by aligning the mRNA sequences of P. canadensis which created by SMRT technology combined with high-throughput sequencing and also considered the Populus trichocarpa genome1. The lncRNAs that have generic exonic overlap with a known transcript was called sense lncRNAs. Contained in the intergenic lncRNA was called long intervening/intergenic non-coding RNAs (lincRNAs). The lncRNAs that have exonic overlap with known transcripts, but on the opposite strand was called antisense lncRNAs, and fall entirely within an intron of a known gene was called sense intronic lncRNAs (Mercer et al., 2009; Ponting et al., 2009; Roberts et al., 2011; Shuai et al., 2014; Wright, 2014; Chen et al., 2015, 2016). Prediction of the target genes of the lncRNAs were cis and trans prediction. To predict the target genes of the lncRNAs, we used BLAST (version 2.9.0+) to search our transcriptome libraries, and select a target sequence that was complementary to the lncRNA, setting E-value = 1e-5 and identity = 80%. Then, the cis-target gene were obtained based on the location relationship between lncRNAs and target genes, and potential target genes (correlation >0.8 or correlation <−0.8, p < 0.05) were selected based on the correlation coefficients of lncRNA and mRNA expression where the target gene was called trans-target gene (Chen et al., 2015; Li et al., 2017; Mao et al., 2017; Lin et al., 2018).
Data Analysis
We used RSEM (Dewey and Bo, 2011) to determine transcript levels from the short-read data. The Illumina data was mapped to PacBio data using LoRDEC software to obtain clean reads. To identify gene expression levels in response to heat stress, the redundant sequences were removed using CD-HIT, and the resulting full-length transcripts were used as a reference sequence (ref) for each gene, and then the clean reads of each sample obtained by short-read were aligned to the ref sequence. Next, we obtained the readcount of each gene by mapping the short-read data to the ref. Further, the readcounts were converted into fragments per kilobase of exon model per million mapped reads (FPKM) values.
Differentially expressed genes (DEGs) were selected using log2 FC≥ 1 or ≤−1 and FDR <0.05 (q-value). Genes with FPKM >0.3 in samples from the two groups were selected for further analysis (Anders and Huber, 2010). All DEGs were mapped to individual terms in the gene ontology (GO) database2 and the number of genes per term was calculated. We then used GOseq software for GO enrichment analysis to identify significant enrichment terms in the DEGs. Analysis of gene regulatory pathways was conducted using the KEGG Pathway database3.
Results
High Quality Full-Length Sequences Were Sufficient for Further Analysis
To reveal the genes and regulatory lncRNAs that respond to heat in P. canadensis, we sequenced and analyzed the transcriptomes of the two RNA pools using the PacBio Sequel platform to obtain accurate long reads. A total of 140,046 polymerase reads and 4,297,179 subreads were produced with SMRT, with an average read length of 2,057 bp (Table 1). To provide more accurate sequence information, we obtained 121,801 circular consensus sequences (CCS) from reads that made at least two full passes through the insert (average read length of 2,490 bp) (Table 1).
Table 1. Summary of reads from third-generation long-read sequencing.
SMRTlink found 114,373 full-length sequences and 111,807 full-length non-chimeric (FLNC) reads (average read length of 2,314 bp) (Table 1). The similar FLNC reads were clustered using the ICE algorithm. The consensus sequence was obtained, and the non-full-length sequence was corrected by arrow software to obtain the consensus sequence, and finally 66,657 polished consensus sequences were obtained (Table 1), with mean length of 2,339 bp.
Then, the cDNA libraries were sequenced following the Illumina HiSeq 2500 platform paired-end protocol. The RNA-seq generated 187,706,638 and 202,259,974 raw data from the pcq and pch samples, respectively, with 183,393,014 (pcq) and 198,641,402 (pch) cleaned reads remaining after trimming (Table 2). The raw reads of Illumina sequencing data were used to correct the SMRT data. After removing redundant and similar sequences, we obtained 86,014,859 nucleotides and 66,657 transcripts with mean length of 2,336 bp. Transcript length distribution results showed that ~1.72, 8.87, 23.6, 42.83, and 22.98% of the transcripts from the corrected isoforms were <500 bp, 500–1,000 bp, 1,000–2,000 bp, 2,000–3,000 bp, and >3,000 bp, respectively. A total of 66,657 corrected consensus sequences were de-redundified using CD-HIT, and 33,840 full-length transcripts were obtained (Table 1).
Table 2. Summary of Illumina sequencing.
High-Quality Full-Length Sequences Improved the Accuracy of LncRNA Prediction
After the corrected consensus sequences were de-redundified using CD-HIT, the 33,840 full-length transcripts were used as reference sequences for the genes. Then the clean reads of each sample obtained by Illumina sequencing were aligned to the ref, and the readcounts of genes were obtained based on the mapping results. The number of mapped reads in the six libraries ranged from 46,427,856 to 61,084,944, and the mapping rates ranged from 86.52 to 89.44% (Supplementary Table S1). In this process, we used RSEM software, and set the parameters of the comparison software bowtie2 in RSEM to end-to-end and sensitive modes, and the others were the default parameters.
For the identification of lncRNAs, we combined the two transcript datasets and selected for transcripts longer than 200 bp and lacking potential coding capability; this predicted lncRNAs was 753 (Figure 1A). We classified these lncRNAs into four groups: 658 sense lncRNAs (87.38%), 41 long intervening/intergenic non-coding RNAs (lincRNAs) (5.44%), 12 antisense lncRNAs (1.59%) and 42 sense intronic lncRNAs (5.58%) (Figure 1B). Analysis of mRNA and lncRNA sequence lengths showed that the average length was 2,564 and 1,568 bp. (Figure 1C). In addition, all 33,840 transcripts (corrected isoforms) were functionally annotated by searching NR, Swissprot, GO, COG, KOG, Pfam, and KEGG databases and a total of 33,676 transcripts (99.5%) were annotated (Figure 1D). We identified the largest number of homologs in P. trichocarpa by comparing the transcript sequences of P. canadensis to the NR database (Supplementary Table S2). Lastly, the Illumina sequencing data were remapped with the SMRT data. The resulting full-length transcripts were used as a ref, and all valid reads were obtained from mapping short-read sequences to ref using RSEM software. After mapping, we converted the readcounts into FPKM to measure the transcript level for each unigene.
Figure 1. Identification of lncRNAs. (A) Venn diagram of lncRNAs predicted by CNCI, CPC, Pfam, and PLEK methods. (B) Proportions of four types of lncRNA. (C) Density and length distributions of mRNAs and lncRNAs in P. canadensis. (D) Seven database annotation result statistics.
Analysis of Differentially Expressed Genes and LncRNAs
To investigate gene expression patterns in P. canadensis, we used FPKM values to measure the transcript level of every unigene and compared the control and the heat-stress-treated groups. This analysis detected 3,493 genes as differentially expressed based on SRMT sequencing (Figure 2A). Among them, 2,601 and 892 were up-regulated and down-regulated, respectively, in the pch group compared with the pcq group. GO enrichment analysis of significantly DEGs showed that the up-regulated genes were associated with 171 biological processes such as response to stimulus, 70 molecular functions such as unfolded protein binding and heat shock protein binding, and 51 cellular components such as transmembrane transporter complex (Figure 3A). Moreover, down-regulated genes were associated with seven biological processes such as response to endogenous stimulus and response to hormones, and six molecular functions such as heme binding, but the cellular components were not found to be significantly enriched (Figure 3B). KEGG pathway enrichment analysis applied on the differentially expressed mRNAs showed that the up-regulated genes were associated with 10 pathways such as protein processing in the endoplasmic reticulum (Figure 4A) and down-regulated genes were associated with 18 pathways such as plant hormone signal transduction (Figure 4B).
Figure 2. Expression profiles of genes and lncRNAs. Hierarchical clustering of all DEGs (A) and lncRNAs (B) in the P. canadensis pcq and pch groups.
Figure 3. GO analysis of the biological function of genes, lncRNAs, and target genes of lncRNAs. GO terms of up- (A) and down-regulated genes (B), up- (C), and down-regulated lncRNAs (D). And up- (E) and down-regulated target genes of lncRNAs (F).
Figure 4. KEGG pathway enrichment analysis of genes, lncRNAs, and target genes of lncRNAs. KEGG pathway of up- (A) and down-regulated genes (B), up- (C) and down-regulated lncRNAs (D). and up- (E) and down-regulated target genes of lncRNAs (F).
Moreover 78 significantly differentially expressed lncRNAs were obtained, including 67 up-regulated lncRNAs and 11 down-regulated lncRNAs (Figure 2B). We identified 923 candidate target genes of lncRNAs, and 72 significantly differentially expressed target genes. GO enrichment analysis of target genes of lncRNAs showed that the up-regulated target genes were associated with 12 biological processes such as protein folding and response to stress, and five molecular functions such as unfolded protein binding and heat shock protein binding, but the cellular components were not found to be significantly enriched (Figure 3C). Moreover, down-regulated target genes were associated with 45 biological processes such as cellular response to stress, 33 molecular functions such as ADP binding, and 26 cellular components such as cytoskeletal part (Figure 3D). KEGG pathway enrichment analysis applied to the target genes showed that the up-regulated target genes were associated with six pathways such as protein processing in the endoplasmic reticulum (Figure 4C) and down-regulated target genes were associated with four pathways such as biotin metabolism (Figure 4D).
LncRNAs Regulate Genes in Response to Heat Stress in P. canadensis
The correlation between lncRNAs and gene expression was calculated by Pearson correlation coefficient, and the interaction network between lncRNAs and target genes was constructed (Figure 5). GO enrichment analysis showed that the up-regulated target genes were associated with five biological processes such as protein folding and response to stress, and one molecular function (unfolded protein binding), but the cellular components were not found to be significantly enriched (Figure 3E). Moreover, down-regulated genes were associated with 33 biological processes such as cellular response to stress, 31 molecular functions such as ADP binding, and 17 cellular components such as mitotic spindle (Figure 3F). KEGG pathway enrichment analysis applied on the differentially expressed target genes showed that the up-regulated genes were associated with three pathways such as protein processing in the endoplasmic reticulum (Figure 4E) and down-regulated genes were associated with one pathway, which was plant-pathogen interaction (Figure 4F). Furthermore, 19 genes in the heat stress-related GO term were annotated, including HSP82, HSP83, DRL27, DNJH, and two predicted protein-coding genes (POPTR_0001s43800g and POPTR_0002s25730g), and these genes are the target gene of four sense lncRNAs (lncRNAPc3, lncRNAPc5, lncRNAPc12, and lncRNAPc14) (Table 3). In addition, four lncRNAs were aligned to the P. trichocarpa transcriptome data by sequence alignment, and the results showed that the genes with the highest ratio of alignment were DRL, UVR8, NDUA9, and HSP83.
Table 3. The DEGs in GO term related to heat stress.
Figure 5. The lncRNA–gene interaction networks. The blue triangle is the candidate lncRNA and the yellow circle is the genes.
Discussion
Full-Length Sequences Identified by SMRT Sequencing in P. canadensis Improved Studies of the Heat Stress Response
The emergence of SMRT sequencing has greatly facilitated the de novo assembly of transcriptomes in higher organisms (Au et al., 2013; Sharon et al., 2013). In Populus tomentosa, short-read sequences obtained from Illumina sequencing showed that the average length of transcripts was 690, 703, and 686 bp in normal, tension, and opposite wood, respectively (Chen et al., 2015), and the average length of transcripts was 671 and 356 bp in Populus euphratica and P. trichocarpa × Populus deltoides, respectively (Qiu et al., 2011; Benjamin et al., 2012). In our study, we first used high-throughput sequencing to obtain high-quality short sequence reads of the P. canadensis transcriptome. Then, we used SMRT sequencing to obtain full-length sequence reads with an average length of 2,339 bp. The sequences obtained by third-generation sequencing were longer than those obtained by second-generation sequencing. In part of lncRNA identification, most of the identified lncRNAs are <2,000 bp in length in our study. The length of lncRNAs were less than the average length obtained by third-generation sequencing, indicating that our method can obtain more accurate lncRNA sequences compared to RNA-seq. This is consistent with the research reported by Wang B. et al. (2016). Furthermore, we explored transcriptomic changes in P. canadensis in response to heat stress by combined SMRT sequencing with Illumina sequencing and demonstrated that this approach provides higher quality reads and more complete assembly compared with using short reads alone (Au et al., 2013; Sharon et al., 2013; Huddleston et al., 2014). We obtained 66,657 complete transcript sequences for P. canadensis by using short reads to correct long reads of SMRT technology sequencing, and greatly improved the accuracy and depth of the study. This is the first study using PacBio Iso-Seq to examine the full-length transcriptome of P. canadensis. The obtained transcriptome may help to further explore P. canadensis genetics and the association of gene expression with phenotypes.
Gene Expression of P. canadensis in Response to Heat Stress
Plants have complex regulatory networks that help them to acclimatize to heat stress (Reddy et al., 2016). Analysis of heat-stress-responsive genes can improve our understanding of the response to heat stress. In this study, we identified 3,493 DEGs that reflect changes in gene expression in P. canadensis heat stress responses and these DEGs provide a foundation for further research. Our GO and KEGG enrichment analysis of these DEGs identified biological processes, metabolic pathways, and biochemical activities that are up- or down-regulated in response to heat, thus providing an in-depth analysis of the transcriptome changes that occur during the heat-stress response of P. canadensis. DEGs associated with high temperature stress in P. canadensis were enriched in the GO terms single-organism metabolic process, oxidation-reduction process, and transmembrane transport biological process. Further study of these DEGs will help us to elucidate the mechanisms involved in heat tolerance and to identify new, potentially useful heat-stress-related genes specific to P. canadensis.
These results were supported by the reports of previous studies (Liu, 2014). Furthermore, we annotated genes that were significantly associated with heat stress with GO terms, and the results showed that many transcription factor-related genes, FtsH protease family genes, HSP family genes, and endopeptidase Clp family genes (including TFB2, FTSH6, HSP90, HSP83, HSP82, and CLPB1) were enriched. Among them, HSP90, HSP83, and HSP82 have been reported to participate in the heat stress response in previous studies (Borkovich et al., 1989; Gkouvitsas et al., 2009; Zhang et al., 2013).
LncRNAs Regulate Genes in Response to Heat Stress in P. canadensis
There are many molecular mechanisms regulating gene expression, including cis-regulatory elements and trans-acting factors such as lncRNAs (Lasky et al., 2014; Albert and Kruglyak, 2015). Our study identified 753 putative lncRNAs in the pch and pcq samples by full-length transcriptome sequencing and classified them into four classes to explore their possible functions. Since the whole genome of P. canadensis has not been sequenced, the genome of P. trichocarpa was used as a reference in our study. The lncRNAs we identified were shorter than the mRNAs, in agreement with previous studies (Pauli et al., 2012). LncRNAs can be significantly differentially expressed in plants under heat stress, for example, Di et al. (2014) found that 245 poly (A)+ and 58 poly(A)– lncRNA have shorter lengths and lower expression levels than mRNAs under heat stress. In our study, 78 lncRNAs showed significant differential expression after heat stress, suggesting that lncRNAs may be involved in the plant high temperature stress responses. However, the regulatory role of these lncRNAs is unclear, and further analysis will be needed to explore this.
LncRNAs can regulate gene expression (Liu et al., 2012). In our study, 923 target genes were obtained by trans prediction. The correlation between lncRNA and gene expression was calculated by Pearson correlation coefficient, and the interaction network between lncRNAs and target genes was constructed. The results of GO enrichment and KEGG pathway analysis indicated that the up-regulated lncRNA target genes were significantly associated with heat stress, suggesting that lncRNAs may regulate target genes in response to heat stress. Further, 19 genes in the heat stress-related GO term were annotated, including HSP82, HSP83, DRL27, DNJH, and two predicted protein-coding genes (POPTR_0001s43800g and POPTR_0002s25730g), and these genes are the target genes of four sense lncRNAs. Among them, HSP82 and HSP83 are the target genes of the sense lncRNAPc5. Further, lncRNAPc5 aligned to P. trichocarpa HSP83, which shows that these lncRNAs identified from full-length transcriptome sequencing may regulate HSP genes by co-expression. These results suggested that lncRNAPc5 may cis-regulate HSP genes via co-expression. Exploration of the functions of the identified lncRNAs will provide further insight into the heat-stress response of P. canadensis.
Data Availability Statement
The raw sequence data reported in this article have been deposited in the Genome Sequence Archive (Genomics, Proteomics and Bioinformatics, 2017) in BIG Data Center (Nucleic Acids Research, 2019), Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under accession numbers CRA002150 and CRA002154 that are publicly accessible at https://bigd.big.ac.cn/gsa
Author Contributions
JC designed the research. JX, MF, ZL, MZ, XL, and YP performed the research. JC analyzed and interpreted the data. JC and JX wrote the paper. YW provided critical revision of the paper. All authors commented on the manuscript.
Funding
This work was supported by Hainan Provincial Natural Science Foundation of China (No. 318MS009), the Opening Project of State Key Laboratory of Tree Genetics and Breeding (Northeast Forestry University, K2018201), National Forest Genetic Resources Platform (No. HD-zkpt-2018001), and the Scientific Research Fund Project of Hainan University (KYQD(ZR)1830).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00249/full#supplementary-material
Footnotes
- ^ http://www.phytozome.net/poplar
- ^ http://www.geneontology.org/
- ^ http://www.genome.jp/kegg/pathway.html
References
Abdel-Ghany, S. E., Hamilton, M., Jacobi, J. L., Ngam, P., Devitt, N., Schilkey, F., et al. (2016). A survey of the sorghum transcriptome using single-molecule long reads. Nat. Commun. 7:11706. doi: 10.1038/ncomms11706
Abdelrahman, M., Elsayed, M., Jogaiah, S., Burritt, D. J., and Tran, L. P. (2017). The “stay-green” trait and phytohormone signaling networks in plants under heat stress. Plant Cell Rep. 36, 11009–1025. doi: 10.1007/s00299-017-2119-y
Ada, N., Norberto, E. P., and Sergio, M. (1987). Early and late heat shock proteins in wheats and other cereal species. Plant Physiol. 84, 1378–1384.
Aimin, L., Junying, Z., and Zhongyin, Z. (2014). A tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics 15:311. doi: 10.1186/1471-2105-15-311
Albert, F. W., and Kruglyak, L. (2015). The role of regulatory variation in complex traits and disease. Nat. Rev. Genet. 16, 197–212. doi: 10.1038/nrg3891
Altschul, S. F., Madden, T. L., Schäffer, A. A., Zhang, J., Zhang, Z., and Miller, W. (1997). Gapped blast and psi-blast: a new generation of protein detabase search programs. Nucleic Acids Res. 25, 3389–3402. doi: 10.1093/nar/25.17.3389
Amor, B. B., Wirth, S., Merchan, F., Laporte, P., Aubenton-Carafa, Y., Hirsch, J., et al. (2009). Novel long non-protein coding RNAs involved in Arabidopsis differentiation and stress responses. Genome Res. 19, 57–69. doi: 10.1101/gr.080275.108.1
Amos, B., and Rolf, A. (2000). The SWISS-PROT protein sequence database and its supplement TrEMBL in 2000. Nucleic Acids Res. 28, 45–48. doi: 10.1093/nar/28.1.45
Anders, S., and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 11:R106. doi: 10.1186/gb-2010-11-10-r106
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: tool for the unification of biology. Nat. Genet. 22, 25–29. doi: 10.1016/S0315-5463(89)70537-X
Au, K. F., Sebastiano, V., Afshar, P. T., Durruthy, J. D., Lee, L., Williams, B. A., et al. (2013). Characterization of the human esc transcriptome by hybrid sequencing. Proc. Natl. Acad. Sci. U.S.A. 110, E4821–E4830. doi: 10.1073/pnas.1320101110
Au, K. F., Underwood, J. G., Lee, L., Wong, W. H., and Xing, Y. (2012). Improving pacbio long read accuracy by short read alignment. PLoS ONE 7:e46679. doi: 10.1371/journal.pone.0046679
Benjamin, P., Emmanuelle, M., Emilie, T., Stéphane, H., Corinne, D. S., Julie, P., et al. (2012). RNA-seq of early-infected poplar leaves by the rust pathogen melampsora larici-populina uncovers ptsultr3;5, a fungal-induced host sulfate transporter. PLoS One 7:e44408. doi: 10.1371/journal.pone.0044408
Bo, W., Elizabeth, T., Michael, R., Tyson, A. C., Ting, H., Yinping, J., et al. (2016). Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing. Nat. Commun. 7:11708. doi: 10.1038/ncomms11708
Borkovich, K. A., Farrelly, F. W., Finkelstein, D. B., Taulien, J., and Lindquist, S. (1989). Hsp82 is an essential protein that is required in higher concentrations for growth of cells at higher temperatures. Mol. Cell. Biol. 9, 3919–3930. doi: 10.1128/MCB.9.9.3919
Chao, Y., Yuan, J., Li, S., Jia, S., Han, L., and Xu, L. (2018). Analysis of transcripts and splice isoforms in red clover (Trifolium pratense L.) by single-molecule long-read sequencing. BMC Plant Biol. 18:300. doi: 10.1186/s12870-018-1534-8
Chen, J., Chen, B., and Zhang, D. (2015). Transcript profiling of Populus tomentosa genes in normal, tension, and opposite wood by RNA-seq. BMC Genomics 16:164. doi: 10.1186/s12864-015-1390-y
Chen, M., Wang, C., Bao, H., Chen, H., and Wang, Y. (2016). Genome-wide identification and characterization of novel lncRNAs in Populus under nitrogen deficiency. Mol. Genet. Genomics 291, 1663–1680. doi: 10.1007/s00438-016-1210-3
Dewey, C. N., and Bo, L. (2011). Rsem: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics 12:323. doi: 10.1186/1471-2105-12-323
Di, C., Yuan, J., Wu, Y., Li, J., Lin, H., Hu, L., et al. (2014). Characterization of stress-responsive lncRNAs in\r, Arabidopsis thaliana\r, by integrating expression, epigenetic and structural features. Plant J. 80, 848–861. doi: 10.1111/tpj.12679
Eddy, S. R. (1998). Profile hidden markov models. Bioinformatics 14, 755–763. doi: 10.1093/bioinformatics/14.9.755
Finn, R. D., Coggill, P., Eberhardt, R. Y., Eddy, S. R., Mistry, J., Mitchell, A. L., et al. (2016). The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 44, D279–D285. doi: 10.1093/nar/gkv1344
Gkouvitsas, T., Kontogiannatos, D., and Kourti, A. (2009). Expression of the HSP83 gene in response to diapause and thermal stress in the moth sesamia nonagrioides. Insect Mol. Biol. 18, 759–768. doi: 10.1111/j.1365-2583.2009.00922.x
Hackl, T., Hedrich, R., Schultz, J., and Forster, F. (2014). Proovread: large-scale high-accuracy pacbio correction through iterative short read consensus. Bioinformatics 30, 3004–3011. doi: 10.1093/bioinformatics/btu392
Havaux, M., Tardy, F., Ravenel, J., Chanu, D., and Parot, P. (1996). Thylakoid membrane stability to heat stress studied by flash spectroscopic measurements of the electrochromic shift in intact potato leaves: influence of the xanthophyll content. Plant Cell Environ. 19, 1359–1368. doi: 10.1111/j.1365-3040.1996.tb00014.x
Heo, J. B., and Sung, S. (2011). Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science 331, 76–79. doi: 10.1126/science.1197349
Hu, T., Sun, X., Zhang, X., Nevo, E., and Fu, J. (2014). An RNA sequencing transcriptome analysis of the high-temperature stressed tall fescue reveals novel insights into plant thermotolerance. BMC Genomics 15:1147. doi: 10.1186/1471-2164-15-1147
Huddleston, J., Ranade, S., Malig, M., Antonacci, F., Chaisson, M., Hon, L., et al. (2014). Reconstructing complex regions of genomes using long-read sequencing technology. Genome Res. 24, 688–696. doi: 10.1101/gr.168450.113
Kant, P., Gordon, M., Kant, S., Zolla, G., Davydov, O., Heimer, Y. M., et al. (2008). Functional-genomics-based identification of genes that regulate Arabidopsis responses to multiple abiotic stresses. Plant Cell Environ. 31, 697–714. doi: 10.1111/j.1365-3040.2008.01779.x
Kong, L., Zhang, Y., Ye, Z. Q., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi: 10.1093/nar/gkm391
Koonin, E. V., Fedorova, N. D., Jackson, J. D., Aviva, R., Jacobs, A. R., and Krylov, D. M. (2004). A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol. 5:R7. doi: 10.1186/gb-2004-5-2-r7
Kotak, S., Larkindale, J., Lee, U., Von, P., Vierling, E., and Scharf, K. D. (2007). Complexity of the heat stress response in plants. Curr. Opin. Plant Biol. 10, 310–316. doi: 10.1016/j.pbi.2007.04.011
Kotoglou, P., Kalaitzakis, A., Vezyraki, P., Tzavaras, T., Michalis, L. K., Dantzer, F., et al. (2009). HSP70 translocates to the nuclei and nucleoli, binds to xrcc1 and parp-1, and protects hela cells from single-strand dna breaks. Cell Stress Chaperones 14, 391–406. doi: 10.1007/s12192-008-0093-6
Lasky, J. R., Des, M. D. L., Lowry, D. B., Inna, P., Mckay, J. K., Richards, J. H., et al. (2014). Natural variation in abiotic stress responsive gene expression and local adaptation to climate in Arabidopsis thaliana. Mol. Biol. Evol. 31, 2283–2296. doi: 10.1093/molbev/msu170
Leone, A., Perrotta, C., and Maresca, B. (2003). “Plant Tolerance to Heat stress: current strategies and new emergent insights,” in Abiotic Stresses in Plants, eds L. S. di Toppi and B. Pawlik-Skowrońska (Dordrecht: Springer).
Li, W., Jaroszewski, L., and Godzik, A. (2002). Tolerating some redundancy significantly speeds up clustering of large protein databases. Bioinformatics 18, 77–82. doi: 10.1093/bioinformatics/18.1.77
Li, Z., Hu, Q., Zhou, M., Vandenbrink, J., Li, D., Menchyk, N., et al. (2013). Heterologous expression of OsSIZ1, a rice SUMO E3 ligase, enhances broad abiotic stress tolerance in transgenic creeping bentgrass. Plant Biotechnol. J. 11, 432–445. doi: 10.1111/pbi.12030
Li, Z., Ouyang, H., Zheng, M., Cai, B., Han, P., Abdalla, B. A., et al. (2017). Integrated analysis of long non-coding RNAs (lncRNAs) and mRNA expression profiles reveals the potential role of lncRNAs in skeletal muscle development of the chicken. Front. Physiol. 7:687. doi: 10.3389/fphys.2016.00687
Lin, Z., Long, J., Yin, Q., Wang, B., Li, H., Luo, J., et al. (2018). Identification of novel lncRNAs in Eucalyptus grandis. Ind. Crops Prod. 129, 309–317. doi: 10.1016/j.indcrop.2018.12.016
Liu, F. (2014). RNA-seq revealed complex response to heat stress on transcriptomic level in Saccharina japonica (laminariales, phaeophyta). J. Appl. Phycol. 26, 1585–1596. doi: 10.1007/s10811-013-0188-z
Liu, J., Jung, C., Xu, J., Wang, H., Deng, S., Bernad, L., et al. (2012). Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis. Plant Cell 24, 4333–4345. doi: 10.2307/23327514
Mao, Y., Liu, R., Zhou, H., Yin, S., Zhao, Q., Ding, X., et al. (2017). Transcriptome analysis of miRNA-lncRNA-mRNA interactions in the malignant transformation process of gastric cancer initiation. Cancer Gene Ther. 24, 267–275. doi: 10.1038/cgt.2017.14
Mathur, S., Agrawal, D., and Jajoo, A. (2014). Photosynthesis: response to high temperature stress. J. Photochem. Photobiol. Biol. 137, 116–126. doi: 10.1016/j.jphotobiol.2014.01.010
Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10, 155–159. doi: 10.1038/nrg2521
Miller, G., and Mittler, R. (2006). Could heat shock transcription factors function as hydrogen peroxide sensors in plants? Ann. Bot. 98, 279–288.
Minoru, K., Susumu, G., Shuichi, K., Yasushi, O., and Masahiro, H. (2004). The KEGG resource for deciphering the genome. Nucleic Acids Res. 32, D277–D280. doi: 10.1093/nar/gkh063
Mora, C., Caldwell, I. R., Caldwell, J. M., Fisher, M. R., Genco, B. M., and Running, S. W. (2015). Suitable days for plant growth disappear under projected climate change: potential human and biotic vulnerability. PLos Biol. 13:e1002167. doi: 10.1371/journal.pbio.1002167
Pauli, A., Valen, E., Lin, M. F., Garber, M., Vastenhouw, N. L., Levin, J. Z., et al. (2012). Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Res. 22, 577–591. doi: 10.1101/gr.133009.111
Place, R. F., and Noonan, E. J. (2014). Non-coding RNAs turn up the heat: an emerging layer of novel regulators in the mammalian heat shock response. Cell Stress Chaperones 19, 159–172. doi: 10.1007/s12192-013-0456-5
Ponting, C. P., Oliver, P. L., and Reik, W. (2009). Evolution and functions of long noncoding RNAs. Cell 136, 629–641. doi: 10.1016/j.cell.2009.02.006
Qiu, Q., Ma, T., Hu, Q., Liu, B., Wu, Y., Zhou, H., et al. (2011). Genome-scale transcriptome analysis of the desert poplar, Populus euphratica. Tree Physiol. 31, 452–461. doi: 10.1093/treephys/tpr015
Reddy, P. S., Chakradhar, T., Reddy, R. A., Nitnavare, R. B., and Reddy, M. K. (2016). “Role of heat shock proteins in improving heat stress tolerance in crop plants,” in Heat Shock Proteins and Plants, eds A. A. A. Asea, P. Kaur, and S. K. Calderwood (Cham: Springer), 283–307. doi: 10.1007/978-3-319-46340-7_14
Roberts, A., Pimentel, H., Trapnell, C., and Pachter, L. (2011). Identification of novel tran-scripts in annotated genomes using RNA-Seq. Bioinformatics 27, 2325–2329. doi: 10.1093/bioinformatics/btr355
Roberts, R. J., Carneiro, M. O., and Schatz, M. C. (2013). The advantages of smrt sequencing. Genome Biol. 14:405. doi: 10.1186/gb-2013-14-7-405
Shamovsky, I., and Nudler, E. (2009). Isolation and characterization of the heat shock RNA 1. Methods Mol. Biol. 540, 265–279. doi: 10.1007/978-1-59745-558-9_19
Sharon, D., Tilgner, H., Grubert, F., and Snyder, M. (2013). A single-molecule long-read survey of the human transcriptome. Nat. Biotechnol. 31, 1009–1014. doi: 10.1038/nbt.2705
Shuai, P., Liang, D., Tang, S., Zhang, Z., Ye, C. Y., Su, Y., et al. (2014). Genome-wide identification and functional prediction of novel and drought-responsive lincRNAs in Populus trichocarpa. J. Exp. Bot. 65, 4975–4983. doi: 10.1093/jxb/eru256
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41, e166–e166. doi: 10.1093/nar/gkt646
Tatusov, R. L., Fedorova, N. D., Jackson, J. D., Jacobs, A. R., Kiryutin, B., Koonin, E. V., et al. (2003). The COG database: an updated version includes eukaryotes. BMC Bioinformatics 4:14. doi: 10.1186/1471-2105-4-41
Treutlein, B., Gokce, O., Quake, S. R., and Südhof, T. C. (2014). Cartography of neurexin alternative splicing mapped by single-molecule long-read mRNA sequencing. Proc. Natl. Acad. Sci. U.S.A. 111, E1291–E1299. doi: 10.1073/pnas.1403244111
Wahid, A., and Close, T. J. (2007). Expression of dehydrins under heat stress and their relationship with water relations of sugarcane leaves. Biol. Plant. 51, 104–109. doi: 10.1007/s10535-007-0021-0
Wahid, A., Gelani, S., Ashraf, M., and Foolad, M. (2007). Heat tolerance in plants: an overview. Environ. Exp. Bot. 61, 199–223.
Wan, L., Zhou, Q., Wang, Y., Wang, W., Bao, M., and Zhang, J. (2015). Identification of heat-responsive genes in carnation (Dianthus caryophyllus l.) by RNA-seq. Front. Plant Sci. 6:519. doi: 10.3389/fpls.2015.00519
Wang, B., Tseng, E., Regulski, M., Clark, T. A., Hon, T., Jiao, Y., et al. (2016). Unveiling the complexity of the maize transcriptome by single-molecule long-read sequencing. Nat. Commun. 7:11708. doi: 10.1038/ncomms11708
Wang, K., Zhang, X., and Ervin, E. H. (2016). “Small heat shock proteins, a Key Player in Grass Plant Thermotolerance,” in Heat Shock Proteins and Plants, eds A. A. A. Asea, P. Kaur, and S. K. Calderwood (Berlin: Springer).
Wright, M. W. (2014). A short guide to long non-coding RNA gene nomenclature. Hum. Genomics 8:7. doi: 10.1186/1479-7364-8-7
Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). Kobas 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. doi: 10.1093/nar/gkr483
Xu, S., Li, J., Zhang, X., Wei, H., and Cui, L. (2006). Effects of heat acclimation pretreatment on changes of membrane lipid peroxidation, antioxidant metabolites, and ultrastructure of chloroplasts in two cool-season turfgrass species under heat stress. Environ. Exp. Bot. 56, 274–285. doi: 10.1016/j.envexpbot.2005.03.002
Yakovchuk, P., Goodrich, J. A., and Kugel, J. F. (2009). B2 RNA and alu RNA repress transcription by disrupting contacts between rna polymerase ii and promoter DNA within assembled complexes. Proc. Natl. Acad. Sci. U.S.A. 106, 5569–5574. doi: 10.1073/pnas.0810738106
Yakovchuk, P., Goodrich, J. A., and Kugel, J. F. (2011). B2 RNA represses TFIIH phosphorylation of rna polymerase II. Transcription 2, 45–49. doi: 10.4161/trns.2.1.14306
Yuan, J., Li, J., Yang, Y., Tan, C., Zhu, Y., Hu, L., et al. (2018). Stress-responsive regulation of long non-coding RNA polyadenylation in Oryza sativa. Plant J. 93, 814–827. doi: 10.1111/tpj.13804
Zhang, J., Li, J., Liu, B., Zhang, L., Chen, J., and Lu, M. (2013). Genome-wide analysis of the Populus Hsp90 gene family reveals differential expression patterns, localization, and heat stress responses. BMC Genomics 14:532. doi: 10.1186/1471-2164-14-532
Keywords: Populus x canadensis Moench, third-generation sequencing, heat stress, lncRNAs, heat shock protein
Citation: Xu J, Fang M, Li Z, Zhang M, Liu X, Peng Y, Wan Y and Chen J (2020) Third-Generation Sequencing Reveals LncRNA-Regulated HSP Genes in the Populus x canadensis Moench Heat Stress Response. Front. Genet. 11:249. doi: 10.3389/fgene.2020.00249
Received: 30 November 2019; Accepted: 02 March 2020;
Published: 07 May 2020.
Edited by:
Rongling Wu, Pennsylvania State University (PSU), United StatesReviewed by:
Junjie Fu, Chinese Academy of Agricultural Sciences, ChinaMilind B. Ratnaparkhe, ICAR Indian Institute of Soybean Research, India
Copyright © 2020 Xu, Fang, Li, Zhang, Liu, Peng, Wan and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jinhui Chen, jinhuichen@hainanu.edu.cn