High-throughput sequencing and transcript assembly
Two rice cultivars, SYD2, a heat-sensitive rice variety, and ZS97B, a heat-tolerant cultivar, were used. Notably, there was less visible damage to seedlings under 2 days of heat stress, but the leaves shrank profoundly after 6 days of stress treatment, especially for SYD2 (Figure S1). We propose that changes in molecular and physiological levels rather than phenotype were the main HS response in the early stage of heat treatment. Three treatments, namely, the control (SYD-0d and ZS-0d), heat stress for 2 days (SYD-2d and ZS-2d) and 6 days (SYD-6d and ZS-6d), were used for transcript profiling. For each condition, the seedlings of three biological replications for each of two varieties were used to construct a total of 18 strand-specific RNA-Seq libraries. Illumina HiSeq X and 150 bp paired-end (PE) reads were generated. After quality control and removal of adaptors with Trimmomatic-0.39, a total of 1,213,356,478 clean reads consisting of 182 billion nucleotides were obtained with an average GC content of 46% (Table S2).
All RNA-seq datasets were mapped to the genome of Oryza sativa.IRGSP-1.0 using Tophat2 with the default parameters except --library-type as fr-firststrand. In total, approximately 1,108,138,248 clean reads were mapped to the reference genome with an average mapped rate of 91.33% and an average of 86.37% reads properly paired mapped (Table S3). Transcript assembly for each library was initially performed using Cufflinks separately, and an average of 102565 transcripts was generated for each library with a high degree of overlap. Merging of the predicted transcript sets resulted in 132054 unique transcripts. PCA suggested that both varieties and heat treatment had an effect on transcriptome differences between the groups and the relatively good repeatability and reliability except the second replicate of ZS97B under 2 days heat stress, Z2b (Figure S2). The sample of Z2b was removed in the subsequent analyses.
Identification and characterization of lncRNAs in rice seedlings
To identify potential lncRNAs in rice under heat stress, the merged transcripts from 18 libraries were used for lncRNA recognition following the pipeline in Figure 1. After removing the known mRNAs, transcripts with lengths less than 200 nt and expression levels lower than 0.5 in all samples, a total of 25,469 potential lncRNA transcripts were identified. Through protein-coding ability inference, 3461, 10787, 9405 and 8809 transcripts were determined to be lncRNAs by CNCI, CPC, PlncPRO and Pfamscan, respectively. With the aim of acquiring highly reliable lncRNAs, the intersection of lncRNAs out of the four programs remained for further analysis. Finally, 2,743 high-fidelity lncRNAs were obtained (Figure 2A).
According to the location of these lncRNAs relative to protein-coding genes, the 2,743 high-fidelity lncRNAs were classified into 1408 intergenic lncRNAs (lincRNAs) (51.3%), 323 antisense lncRNAs (lncNATs) (11.8%), 960 sense lncRNAs (35.0%), and 52 intronic lncRNAs (1.9%) (Figure 2B). The lncRNAs were distributed across the 12 rice chromosomes relatively evenly (Figure 2C), consistent with those observed in other species, Arabidopsis thaliana (Severing et al. 2018), maize (Li et al. 2014) and Gossypium hirsutum (Deng et al. 2018). The length of the majority of lncRNAs was concentrated at approximately 1000 nt, with an average of 1324 nt and a median length of 969 nt, shorter than that of mRNA 1591 nt and 1419 nt (Figure 2D). With regard to the number of exons, approximately 48.9% of lncRNAs consisted of a single exon. A total of 39.7%, 10.8% and 1.3% of lncRNAs harboured 2-4, 5-9 and over 10 exons, respectively. In contrast, there were fewer overall exons in lncRNAs than in mRNAs. mRNAs with single, 2~4, 5~9 and over 10 exons accounted for 27.0%, 37.6%, 23.3% and 12.2%, respectively (Figure 2E). The overall expression trend was similar across the 6 groups between lncRNAs and mRNAs. However, the expression levels of lncRNAs ranged from 1 to 5 FPKM, which was lower than that of mRNAs (Figure 2F). Similar expression features of lncRNAs vs. mRNAs have been found in Arabidopsis (Liu et al. 2012), cotton (Wang et al. 2015), tomato (Wang et al. 2018) and maize (Li et al. 2014), indicating the conservation of lncRNAs in these characteristics.
Identification and annotation of target-coding genes of lncRNAs in rice seedlings
It has been suggested that lncRNAs can perform their functions through diverse mechanisms. We identified target genes of lncRNAs in accordance with their functional ways, mainly including cis action, ceRNA, antisense transcripts and precursors of microRNA (Chekanova, 2015). Among these lncRNAs, cis-regulation suggests that some lncRNAs can regulate the expression of their proximal protein-coding genes. A total of 5,406 protein-coding genes were detected within 10 kb upstream and downstream of 2,435 lncRNAs, accounting for 88.7% of the total number of lncRNAs, suggesting that most of the lncRNAs were located in coding gene enrichment regions. The steady-state FPKM values of lncRNAs and their adjacent protein-encoding genes were subsequently used for co-expression analysis based on Pearson correlation coefficients (|r| > 0.5 and P < 0.05). Overall, 1278 lncRNA-mRNA co-expression pairs were identified, involving 926 lncRNAs and 1079 protein-coding genes. These genes are potential targets of lncRNAs by cis-action (Table S4).
Some antisense lncRNAs have been found to be able to bind to their corresponding sense mRNAs to affect the posttranscriptional modification or stability of the mRNAs (Li et al. 2015a). The thermodynamic minimum free energy between lncRNA and its sense chain mRNA calculated by LncTar showed 237 lncRNAs complementary to their sense strand mRNA. More than 80% of these lncRNAs had very high levels of interaction with their sense mRNAs (low: ndG<-0.1; high: ndG<-0.15; and very high: ndG<-0.2) (Table S5), implicating potential roles of these lncRNAs in the posttranscriptional regulation of corresponding sense-coding genes.
Competing endogenous RNAs (ceRNAs) can be targeted by specific miRNAs via target mimics due to sharing miRNA recognition elements (MREs) and thereby can regulate the interaction between miRNAs and their target mRNAs. When ceRNA is silenced, mRNA is degraded by the silencing complex mediated by microRNA (RISC). When ceRNA is transcribed, it can compete with the RISC complex to reduce the inhibitory function of microRNA and increase the expression of target mRNAs. Thus, potentially positive regulation exists between ceRNAs and the target genes of their common microRNAs (Thomson & Dinger 2016). We predicted 915 lncRNA-microRNA pairs and 8692 microRNA-gene target pairs by using psRNATarget (Dai et al. 2018). Based on this, 26097 potential ceRNA pairs of lncRNA-microRNA-gene were generated. We further calculated the expression correlation of lncRNAs and coding genes at the steady-state level in the same lncRNA-microRNA-gene suite. Finally, 2134 ceRNA pairs of "lncRNA-microRNA-gene" were obtained with a significant positive expression correlation between lncRNA and gene" (Table S6). The ceRNA network contained 1204 nodes, including 317 lncRNAs, 671 genes and 216 microRNAs (Figure S1A). Multiple lncRNAs and multiple genes can share one or a few microRNAs (Figure S3B, S3C), while a single microRNA can regulate multiple genes and lncRNAs simultaneously (Figure S3D).
Another mode of lncRNA functions is serving as a precursor of microRNAs to be spliced in vivo to form microRNAs to further regulate other coding genes (Wilusz et al. 2009). We compared sequences of lncRNAs identified in the present study and rice microRNAs. Eighty-nine lncRNAs were highly homologous to 121 microRNAs from 36 families. Subsequently, the target protein-coding genes of the 121 microRNAs were predicted. Finally, 89 lncRNAs (predicted precursors of microRNAs), 121 microRNAs and 675 genes formed 6705 lncRNA-microRNA-gene correlations (Table S7).
GO enrichment analysis showed that 789 (40.1%) target genes were annotated in the category of biological processes. Among these genes, 534, 532, 218 and 206 were involved in metabolic processes, cellular processes, the response to stimuli and biological regulation, respectively. A total of 650 (33.0%) target genes were enriched in the category of cell components, among which 464 genes were in the category of cell and cell part, followed by 356 genes enriched in membrane and 341 of those in organelle. A total of 951 (48.3%) target genes were annotated in the category of molecular function. A total of 682 of them were enriched in the category of binding activity and 482 in catalytic activity. In addition, some target genes were classified into the categories of transporting activity, nucleic-acid binding, transcription factor activity and electronic carrier activity (Figure S4A). KEGG pathway annotation showed that the target genes of lncRNAs were involved in metabolic pathways (specifically, carbohydrate metabolism, amino-acid metabolism, energy metabolism and so on), genetic information-processing translation, folding, classification and degradation of proteins and environmental information responses (Figures S4B, S5, S6, S7, S8). These results indicate the extensive pathways in which lncRNAs are involved.
Expression pattern of lncRNAs responding to heat stress
A total of 128 differentially expressed lncRNAs (DELs) at the steady-state level were detected in SYD2, with 73 DELs after 2 days of stress and 93 DELs after 6 days of stress compared with the control. Thirty-eight of them were shared at both time points of stress. A total of 134 DELs were detected in ZS97B, including 94 and 95 at days 2 and 6 compared with the control, respectively. Thirty-eight of them were shared at both stress-time points. Overall, 231 lncRNAs potentially responsive to high-temperature stress were detected, and 31 of them were responsive to heat stress in both varieties. A total of 103 and 97 lncRNAs were responsive to heat stress in ZS97B and SYD2, respectively (Figures 3A and 4A, 4C, 4E). To validate their expression patterns, a few lncRNAs and target-coding genes were selected to measure their expression levels with quantitative real-time reverse transcription-PCR (qRT-PCR). We selected the lncRNA and target gene pairs to cover different predicted regulation patterns and different expression trends. In detail, the Os01g0160800:TCONS_00018499 pair showed the ceRNA relationship mediated by osa-miR1439, and they revealed heat-induced expression in a heat-sensitive variety. osa-miR1439 has been indicated to play a specific role in heat-stress regulation (Mueller et al. 2017). The Os01g0326100:TCONS_00043075 pair also showed the ceRNA relationships mediated by osa-miR1436, but they are both downregulated in each variety. The expression of osa-miR1436 has previously been demonstrated to be induced in response to heat stress in both tolerant and susceptible rice varieties (Mangrauthia et al. 2017b). Os03g0781700:TCONS_00070624, Os05g0582300:TCONS_00088405 and Os06g0195800:TCONS_00100154 represent cis-regulation relationships and exhibit opposite expression trends, similar repressed expression and similar upregulated trends under heat stress, respectively. The relative expression trends for the lncRNAs and coding genes were mostly consistent between these two methods (Figure 5).
Interaction network of heat-responsive lncRNAs with their target genes
The interaction networks between heat stress-responsive lncRNAs (HRLs) and their targets were constructed for three groups that were responsive to heat stress in ZS97B and SYD2 and shared in both. The interaction relationships of lncRNAs and target genes were highlighted as cis-action, antisense, precursor of miRNA and ceRNA. It is noteworthy that, except for the target genes screened through the former four approaches mentioned above, we calculated the expression Pearson correlation coefficients of the steady-state FPKM values between HRLs and the restprotein-encoding genes. Twenty-four lncRNA-gene pairs were discovered for heat-responsive lncRNAs shared in both varieties, and 88 and 86 lncRNA-gene pairs were discovered for lncRNAs responsive to heat stress in the ZS97B and SYD2 varieties (|r| > 0.9 and P < 0.05), respectively. These lncRNA-gene pairs were not detected by the four above target-gene prediction methods, cis-targets, antisense target ceRNA pairs or precursors of microRNAs. We designed them as trans-action pairs and proposed that other regulatory mechanisms exist between them.
For lncRNAs responsive to heat stress shared in both varieties, 238 lncRNA-mRNA pairs were identified, including relationships by 71 ceRNAs, 4 antisenses, 96 cis-actions, 43 precursors of miRNA and 24 trans-actions (7 negative and 17 positive relations) (Figures 6A, 7, 8, S9, Table S8). For lncRNAs responsive to heat stress in ZS97B, 173 lncRNA-mRNA pairs were identified, including relationships by 11 ceRNAs, 70 cis-actions, 4 precursors of miRNAs and 88 trans-actions (21 negative and 67 positive relations) (Figures 6B, S10, Table S9). For lncRNAs responsive to heat stress shared in SYD2, 168 lncRNA-mRNA pairs were identified, including relationships by 9 ceRNAs, 71 cis-actions, 2 precursors of miRNAs and 86 trans-actions (13 negative and 73 positive relations) (Figures 6C, S11, Table S10).
Functional analysis of heat-responsive lncRNAs by their target genes
The functions of lncRNAs can be predicted indirectly by annotation of their target genes. The 380 DEL target genes were subjected to KEGG pathway enrichment analysis, and 178 of them were annotated (Table S11). The significantly enriched KEGG classes refer to amino-acid metabolism, biosynthesis of other secondary metabolites, carbohydrate metabolism, energy metabolism, protein folding, sorting and degradation. The extensive regulatory roles between lncRNAs and protein-encoding genes highlighted the importance of lncRNAs in the plant response to heat stress, mainly by participating in oxidation-reduction process association, photosynthesis, translation processes and protein processing (Figure 3B).
We further performed GO enrichment analysis for HRL target genes in different groups. For lncRNAs responsive to heat in both varieties, GO enrichment analysis was conducted for 187 target genes. GO terms of photosynthesis, energy-coupled proton transport, ATP synthesis, mitochondrial proton transport, and active ion transmembrane transporter activity were significantly enriched. We proposed that these GO terms are general heat-response processes and DELs and that their target genes enriched in these GO terms are important regulons in response to heat stress in both varieties, regardless of their susceptibility to stress (Figure 4B). For lncRNAs responsive to heat in ZS97B, GO enrichment analysis was conducted for a total of 99 target genes. Significantly enriched GO terms included nucleoside triphosphate metabolic processes, dehydrogenase (NAD+), Hsp70 protein binding and so on (Figure 4D). These results indicated that nucleoside triphosphate metabolic process-, dehydrogenase (NAD+)-, and Hsp70-related lncRNAs play vital roles in ZS97 tolerance to heat stress. For lncRNAs responsive to heat in SYD2, GO enrichment analysis was conducted for 110 target genes, and the significantly enriched GO terms included xyloglucan transferase activity, protein dimerization activity, apoplasts, and RNA polymerase II transcription factor complexes (Figure 4F). We suggested that the pathways harbouring these lncRNAs and genes were significantly affected in the sensitive variety SYD2 under heat stress. The functional annotation of the target genes indicated that lncRNAs played important roles in plant-stress responses in multiple ways.
The critical functions of heat-responsive lncRNAs through interaction with coding genes
In the present study, the various interplay modes between HRLs and target genes were shown in the interaction network. The significance of some lncRNAs and coding genes were highlighted as their hub positions in the network. For instance, lncRNA TCONS_00043075 is a target of osa-miR1436 (Figure 6A). TCONS_00043075, along with 17 other protein-coding genes, forms ceRNA pairs mediated by osa-miR1436. The roles of these genes refer to multibiological processes and molecular functions, including response to oxidative stress (Os01g0326100), response to abiotic stimulus (Os12g0264500, enoyl-CoA hydratase/isomerase family protein), mRNA capping (Os02g0780600), protein binding (Os06g0167500, leucine-rich repeat containing), zinc-ion binding (Os06g0340200, zinc finger), NB-ARC domain containing (Os11g0266500) and kinase-like domain (Os11g0513700)-containing transposon proteins (Os12g0292200). The expression of osa-miR1436 has previously been demonstrated to be induced in response to heat stress in both tolerant and susceptible rice varieties (Mangrauthia et al. 2017b). Consistently, TCONS_00043075 was downregulated in both varieties at high temperature.
LncRNA TCONS_00092993 was predicted to be a precursor of microRNA osa-miR1850, which was predicted to target 39 protein-coding genes (Figure 6A, Figure 7). osa-miR1850 was demonstrated to be downregulated in a heat-tolerant variety, N22, shortly after heat-stress treatment (Mangrauthia et al. 2017a). The osa-miR1850 precursor, TCONS_00092993, was induced under heat stress in both varieties in the present study, indicating complex regulatory mechanisms among these RNA molecules. The target genes of osa-miR1850.1 were involved in stress-relevant pathways. For instance, Os01g0337900 was associated with cellular redox homeostasis (GO: 0045454), Os02g0711300 encoded a heat-shock protein, Hsp20, Os03g0594100 was involved in the oxidation-reduction process, Os05g0528900 encoded a ribosomal protein, and Os06g0113950 encoded a RuBisCO large subunit related to photosynthesis. These results suggested that TCONS_00092993 may indirectly regulate the expression of multiple stress-responsive genes by serving as a precursor of osa-miR1850, being spliced and subsequently regulating the expression of target genes in rice under heat stress.
The lncRNA, TCONS_00100154 (Figure 6A), and its adjacent downstream gene (Os06g0195800, heat-shock protein DnaJ) were hardly expressed under optimal conditions in both varieties ZS97B and SYD2 but were induced after heat stress treatment. TCONS_00005591 (Figure S9) and its adjacent genes (Os01g0719100, OsRZFP34) were co-induced by heat stress. It has been reported that OsRZFP34 regulates stomatal openings in rice under heat and drought stresses (Hsu et al. 2014). All of these results indicated that lncRNAs could play vital regulatory roles during the adaptation of rice to heat stress.
In addition to lncRNAs in the hub position of the network, some lncRNAs co-interact with single coding genes, which are at the centre of the network. For example, the expression of Os01g0101600, an RNA-binding motif, negatively correlated with 7 lncRNAs (Figure 6C), whose cis-target genes were Os02g0216600 (Protein ROOT HAIR DEFECTIVE 3), Os12g0424180 (ATP synthase, LOC_Os12g23610.1), Os06g0547100 (peroxidase activity, response to oxidative stress) and Os09g0321900 (acid-amino acid ligase activity). Os01g0100100 (Figure 6B) was positively correlated with 17 lncRNAs, which in turn regulated other coding genes in cis mode or as a precursor of a miRNA. Os11g0576100 (Figure 6B) (LOC_Os11g36760 transmembrane receptor, GO:0006950 response to stress) was the cis-target of four lncRNAs. We proposed that coding genes with regulatory functions such as RNA-binding activity can regulate heat-responsive lncRNAs, which subsequently regulate important stress-responsive downstream genes.
Among the lncRNAs in response to heat stress in SYD2, TCONS_00001878 with Os01g0104900, TCONS_00018499 with Os01g0160800 and TCONS_00030558 with Os01g0196800 form ceRNA pairs, where osa-miR1439 functions as a mediator (Figures 6C, S11). The expression of osa-miR1439 was induced during high temperature, suggesting a specific role of osa-miR1439 in heat-stress regulation (Mueller et al. 2017). Os01g0104900 and Os01g0160800 encode transferase activity and the ribosome inactivating protein, OsRIP1, respectively. The translation inhibitory activity and the activity of RIPs both increased in leaves when subjected to heat or osmotic stress. It has been proposed that a physiological role of RIPs could be to intervene in the death of plant cells under heat or osmotic stress (Stirpe et al. 1996). These findings together suggested that lncRNAs could regulate plant responses to heat stress through a ceRNA mode to participate in lncRNA-osa-miR1439-gene (such as ribosome-inactivating proteins) circuits.
Transcription factors were also involved in the response to heat stress in SYD2. Gene Os01g0104200 encodes a DNA-binding, no apical meristem (NAM), domain-containing protein. Its homologue in Arabidopsis AtRD26, encoding a NAC transcription factor, was annotated RESPONSIVE TO DESICCATION. The Os01g0104200 gene was positively co-expressed with 17 lncRNAs (Figure 6C). We speculate that this NAC TF functions as an upstream regulatory factor regulating these lncRNAs. LncRNAs then regulate other stress-related genes in different fashions. For example, Os07g0628500 was regulated by TCONS_00109028 and TCONS_00124705 targets Os09g0263933 (similar to the flavonoid 3-monooxygenase.oxidation-reduction process), TCONS_00109958 targets Os07g0105600 (Figure 6C: Illustration of the cis-relationship between them) (Photosystem II oxygen-evolving complex protein PsbQ family protein), TCONS_00072956 targets Os04g0271200, TCONS_00072957 targets Os04g0271200, TCONS_00050981 targets Os02g0211200, TCONS_00017573 targets Os10g0124900, TCONS_00119582 targets Os08g0565200, TCONS_00030589 targets Os11g0219000 (Peptidase S16, a N-terminal domain-containing protein), TCONS_00011753 targets Os01g0381325, TCONS_00011754 and TCONS_00011755 targets Os01g0381325 in a cis fashion. TCONS_00039696 serves as a target mimic of osa-miR156k, regulating Os01g0242400, whereas TCONS_00082845 was predicted as a precursor of miRNA osa-miR5505, regulating Os01g0348700 (similar to 60S ribosomal protein L23a), and TCONS_00078463 as a precursor of osa-miR530-3p, regulating Os01g0192600 (Figure 6C).
Colocalization of heat-responsive lncRNAs with rice heat stress-related QTLs
QTLs for seedling growth under heat stress were collected from the study of Kilasi et al. (2018). In this study, QTLs were identified for seedling growth in response to high temperature based on a recombinant inbred line population constructed by crossing heat-tolerant "N22" and heat-susceptible "IR64" varieties. The traits and materials were close to those in this study. The relevant QTL intervals and their chromosome physical positions were retrieved from the study and compared with the positions of heat-responsive lncRNAs identified in the present study. 10 DELs were colocalized with 5 QTLs (Figure 9). No target-coding genes were predicted for lncRNA TCONS_00067817, located within QTL slpc3.1 on chromosome 3, and TCONS_00022346, within slpc10.2 on chromosome 10, where QTL slpc was denoted as "shoot length under heat stress expressed as percentage of control". Seven DELs and their cis-target genes were colocalized within the same QTLs. For example, lncRNA TCONS_00070624 was localized on chromosome 3 QTL rlpc3.1, denoted as "root length under heat stress expressed as percentage of control". Six cis-target-coding genes of TCONS_00070624 were predicted, including Os03g0781100, Os03g0781200, Os03g0781300, Os03g0781400, Os03g0781600 and Os03g0781700. Remarkably, Os03g0781100, encoding a pentatricopeptide repeat (PPR) protein, has been demonstrated to bind RNA molecules and to function in diverse RNA metabolic pathways (Kobayashi et al. 2012). We speculated that Os03g0781100, encoding a PPR protein, might bind TCONS_00070624 to contribute to the response to high temperature. TCONS_00077715 and its cis-target gene, Os04g0120350, were mapped to QTLs rlpc4.1 and slht4.1 (shoot length under heat stress). Os04g0120350 encodes a protein similar to H0716A07.11 involved in the biological process of proteolysis (GO:0006508), indicating that TCONS_00077715 and Os04g0120350 cooperate to function in proteolysis induced by high temperature. TCONS_00084651 and Os05g0194900 were positioned within QTL slpc5.1 on chr5. Within QTL slpc10.2 on chromosome 10, 6 heat-responsive lncRNAs were located, among which TCONS_00018499 was adjacent to the cis-target-coding gene, Os10g0329400, and TCONS_00018800 to Os10g0368902, TCONS_00018917 to Os10g0369000 and TCONS_00018921 to Os10g0387000. Interestingly, Os10g0368902 also encodes a PPR protein, suggesting the vital roles of lncRNA- and PPR-coding gene interactions in plants coping with heat stress. Additionally, 3 DELs within QTLs had some target genes by the mode of ceRNA. Both TCONS_00018455 and Os01g0720600 were targeted by microRNA osa-miR5071. TCONS_00077715 together with Os05g0248301 and Os11g0227600 form the ceRNA pairs through osa-miR2090. TCONS_00018499 can serve as a sponge RNA of osa-miR1439 for 7 protein-coding genes (Figure 9). Notably, some of them encode heat response-related proteins; for instance, Os01g0160800 is relevant to protein synthesis, and Os07g0601900 encodes a protein similar to NADPH HC toxin reductase.
Enriched sequence motifs within heat-responsive lncRNAs
RNAs bound and regulated by proteins with similar functions usually contain conserved sequence motifs (Ray et al. 2013). LncRNAs induced under the same stress may be regulated by similar transcription factors or other regulated proteins. It has previously been demonstrated that stress-related sequence motifs were enriched in stress-responsive lncRNAs (Di et al. 2014; Zhang et al. 2019). In this study, we predicted the sequence motifs of lncRNAs responsive to increasing temperature. For the 231 heat-responsive lncRNAs, MEME was used for conserved sequence motif searching with the parameters mentioned in the Methods section. Ten motifs were found in all, with the most site occurrences of 125 and the least site occurrences of 6. Five motifs with the top occurrence number and lowest E-values were predominantly enriched with CGG, CUC, CG, AAA and UUU (Figure 10). To further verify motif enrichment in HRL body sequences, three background datasets, namely, HRL promoters, non-HRL body sequences and non-HRL promoters, were selected to compare the occurrence of motifs. By searching the motifs with MAST software in background sequences and Fisher's exact test, we found significant enrichment differences of motifs 2, 4, and 5 between HRL bodies and non-HRL bodies; motif 1 between HRL bodies and HRL promoters; and motifs 1, 2, 3, and 4 between HRL promoters and non-HRL promoters (Figure 10B). These results suggested the importance of lncRNAs with these motifs in the response to heat stress.