The Signicance of Heat Responsive Long Noncoding RNAs in Rice Discovered by Integrating Expression, Network and Genetic Analyses

Long noncoding RNAs (lncRNAs) play vital roles in plant responses to environments. However, the systematic identication of lncRNAs in rice under heat stress remains to be fully performed. We performed steady-state strand-specic RNA-seq for rice seedlings under heat stress to predict lncRNAs. In total, 2743 lncRNAs were predicted, and their expression proles in response to heat treatments were provided. We identied 231 differentially expressed lncRNAs (DELs) in their steady-state level under heat stress, including 31 DELs common to both varieties and 103 and 97 specic to ZS97B and SYD2, respectively, all dened as heat-responsive lncRNAs (HRLs). The target-coding genes of HRLs were predicted through cis-action, ceRNAs, precursors of miRNA, antisense modes and other trans-action targets. GO and KEGG enrichment analyses of HRL targets revealed functions in which HRLs were involved. The interaction network between HRLs, target genes and relevant miRNAs was constructed. The HRLs and their targets were integrated with publicly available QTLs for rice seedling growth under heat stimulus, and 10 HRLs and 12 target genes were co-localized with 5 heat-relevant QTLs. Sequence motif analysis found a few signicant motifs enriched in 231 HRL sequences, and the differential signicance of enriched motifs between HRLs and background sequences was tested. Our ndings provide new insights into the characterization of lncRNAs in terms of the heat response and resources for further investigations of plant heat tolerance improvement. QTLs for seedling relevant QTL intervals and their chromosome physical positions were retrieved from the study and compared with the positions of heat-responsive lncRNAs. Those QTLs and lncRNAs with overlapped positions were selected and the QTLs, overlapped lncRNAs and their targets were then mapped on the chromosome with MapChart 2.3 (Voorrips 2002). strategies Elucidating the roles of lncRNAs in regulating rice adaptability to heat stress in study will provide a novel molecular basis to cope with high temperature. Using steady-state strand-specic RNA-seq, we identied 231 heat-responsive lncRNAs in rice seedlings. An interaction network of HRLs and target-coding genes was constructed to predict the potential roles of HRLs in heat-stress responses. GO enrichment analysis of HRL targets highlighted the general role of lncRNAs in photosynthesis, ATP synthesis, mitochondrial-energy metabolism, transmembrane transporter activity and heat-shock protein Hsp20 at high temperature, independent of variety. We also identied the tolerant variety of ZS97B-specic responsive lncRNAs, which were predicted to potentially interact with coding genes in the nucleoside triphosphate 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 were detected. The results indicated the importance of these motifs in specic lncRNAs functioning in response to heat stress. However, further testing in a greater number of diverse datasets and experimental validation are required to strengthen the ndings here. We anticipate further sequence motif elucidation of abiotic stress-responsive lncRNAs to provide a framework for future biotechnological applications of lncRNAs in crop improvement in abiotic stress tolerance, such as target modication relevant to important motifs.


Introduction
Global climate change poses a great challenge to plant productivity and grain security (Lancaster & Humphreys 2020). Increasing temperature is one of the major adverse environmental factors that in uences plant growth, development and productivity through negative effects on the morphological, physiological, metabolic and molecular levels in plants (Wilkins et al. 2016). Generally, the negative effects of heat stress on plants include impairing photosynthetic activity, producing reactive oxygen species (ROS), decreasing cell expansion and division, membrane uidity, reducing the photosynthesis rate, DNA breakage, protein denaturation and so on (  . The antisense lncRNA COOLAR can physically associate with Arabidopsis thaliana Flowering Locus C (FLC) and accelerate transcriptional shutdown of FLC during vernalization through switching of chromatin states (Csorba et al. 2014). LncRNAs can competitively bind to nuclear speckle, RNA-binding protein (NSR) as competitor long noncoding RNAs to disturb the splicing of NSR mRNA targets (Bardou et al. 2014). A SNP in the promoter region of long-day-speci c, male fertility-associated RNA (LDMAR) can increase methylation of the putative promoter region of LDMAR, which reduces the transcription of LDMAR, resulting in premature programmed cell death (PCD) in developing anthers (Ding et al. 2012). An Arabidopsis lncRNA, DRIR, has been elucidated as a positive regulator of the plant response to drought and salt stress (Qin et al. 2017). LncRNA MuLnc1 has been demonstrated to participate in the response to salt and drought stresses through the interaction with miR3954 and Calmodulin-like protein gene CML27 in mulberry (Gai et al. 2018).
Along with the functional characterization of several lncRNAs in plant stress-responsive pathways, a comprehensive survey of lncRNAs responsive to diverse environmental stresses has been conducted in a variety of species. Genome-wide analysis of natural antisense transcripts (NATs) in Arabidopsis identi ed 1392 spatial-and developmental-speci c NAT pairs in response to light and indicated a strong correlation between light-regulated NAT pairs and histone acetylation on their promoter and gene-body regions (Wang et al. 2014a). LncRNA-mRNA co-expression pro les in response to drought stress have been explored in a range of plant species, including wheat (Cagirici et al. 2017), foxtail millet (Qi et al. 2013), Populus trichocarpa (Qi et al. 2013; Shuai et al. 2014) and rice (Singh et al. 2017). Identi cation of lncRNAs responsive to salt stress has been performed in Gossypium hirsutum (Deng et al. 2018). As an important environmental factor impairing plant growth and productivity, the effect of high temperature on lncRNA expression has been pro led in various species, such as Arabidopsis (Severing et al. 2018), Chinese cabbage (Song et al. 2016) and wild banana (Liu et al. 2018a). However, how lncRNAs respond to heat stress and regulate target gene expression in rice remains to be elucidated.
In this study, we systematically identi ed and characterized lncRNAs in two rice varieties with contrasting heat tolerance with the strand-speci c RNA-seq (ssRNA-seq) technique. Through analysing the expression pro le of lncRNAs and their interaction patterns with protein-coding genes and small RNAs, the sequence motif within HRL sequences, and the overlap of HRLs and heat adaptation-related QTLs, we aimed to draw a holistic picture of how lncRNAs respond to heat stress in rice seedlings and screen out the signi cant lncRNAs and their key sequence motifs to promote the breeding of rice tolerant to heat stress.

Methods
Plant materials, growth conditions and stress treatment Initially, two rice cultivars, SYD2 (O. sativa L. ssp. japonica) and ZS97B (Oryza sativa L. ssp. indica) were used. SYD2 is a heat-sensitive rice variety, whereas ZS97B is a heat-tolerant cultivar, as the maintainer line for a number of elite hybrids widely cultivated in China. Germinated seeds were sown in PVC pots (9 cm diameter and depth) lled with soil under 80% relative humidity in a growth chamber with a 14 h/25°C day and 10 h/15 °C night cycle. Germinating seedlings were kept in these optimal conditions for two weeks. Heat stress was imposed by changing the growth chamber temperature with 14 h/42 °C day and 10 h/20 °C night cycles. On the 2nd day and 6th day of heat-stress treatment, the leaf tissue of 5 plants from each pot was harvested and pooled together for each cultivar as the short and long stress treatments. Leaf tissue under optimal growth conditions was harvested as a control. All harvested leaves were immediately stored at -80 °C for further analysis. Three biological replicates were applied for each condition and cultivar.
Construction of the lncRNA library and Ribo-minus sequencing . Brie y, coding potential scores smaller than 0 were considered noncoding RNAs using CPC2 with default parameters. CNCIs were used with the parameters "-m: pl" and "score < 0". PlncPRO was used to retain the transcripts with "label=0". The aminoacid sequences of the intersecting transcripts ltered by CPC2, CNCI and PlncPRO were compared with the pfamA database using PfamScan, and transcripts with E-values smaller than 0.0001 were ltered out. The rest of the transcripts were considered candidate lncRNAs. Based on the location of lncRNAs relative to protein-coding RNAs, lncRNAs were classi ed as intergenetic lncRNAs (lincRNAs), intronic lncRNAs, sense lncRNAs, and antisense lncRNAs, which were distinguished by Cuffcompare as class codes of u, i, o and x, respectively (Chekanova 2015). The number, length and expression of lncRNA exons were counted by an in-house Python script, and visualization was carried out in R packages. In this study, all identi ed lncRNAs were compared with known rice lncRNAs in the lncRNA databases GREENC (Paytuvi Gallart et al. 2016) and CANTATAdb (Szczesniak et al. 2016), and lncRNAs with an E-value < 1e-10 and sequence coverage > 55% were de ned as known lncRNAs.
Gene expression quali cation and differential expression analysis

Prediction of lncRNA targets
LncRNAs are implicated in exerting critical functions by interacting with protein-coding gene expression in different ways. First, plant antisense lncRNAs can play their roles via interactions with sense mRNAs. To evaluate the interaction between antisense lncRNAs and their sense mRNAs, we used LncTar to predict their complementary binding ability, with ndG<-0.15 as the threshold (Li et al. 2015a). Second, lncRNAs were reported to be able to regulate the transcription of their neighbouring genes by cis-action. To this end, we searched for protein-coding genes located in the 10 kb upstream and downstream regions of lincRNAs. Furthermore, we calculated the co-expression correlation between lncRNAs and adjacent genes based on the Pearson correlation coe cient. Those adjacent genes were co-expressed with their neighbouring lncRNAs with |R| > 0. 5 and P < 0.05 were considered cis-target genes of lncRNAs (Luo et al. 2016). Third, lncRNAs can function in biological processes as precursors of microRNAs (Shumayla et al. 2017a). We compared lncRNAs identi ed in this study with rice microRNAs in the Microbase database (http://www.mirbase.org). lncRNAs with coverage of known microRNAs greater than 90% and Evalue < 1E-10 were identi ed as precursors of microRNAs. Fourth, some lncRNAs can function as competitive endogenous RNAs (ceRNAs) (Cesana et al. 2011;Xu et al. 2016). psRNA Target was used to predict the targets of microRNAs (lncRNAs and protein-coding genes) (Dai et al. 2018), and the ceRNA correlation was determined if the "lncRNA-mRNA" pairs were co-expressed with correlation R > 0.5 and P < 0.05 and were targeted by the same miRNA via target mimicry ).

Analysis of DEGs and lncRNA function
To inspect the functions of DEGs, the RAP-DB database (https://rapdb.dna.affrc.go.jp/) was used to convert rice gene IDs from RAP IDs to MSU IDs. Different genes were annotated and classi ed based on the RGAP (http://rice.plantbiology.msu.edu/index.shtml) database. The PlantDB (http://planttfdb.cbi.pku.edu.cn/) database was used to identify members of the transcription factor family in rice (Jin et al. 2017). Gene ontology (GO)-based enrichment tests were conducted using Singular Enrichment Analysis (SEA) in the AgriGO toolkit (http://bioinfo.cau.edu.cn/agriGO/analysis.php). Oryza sativa was selected as the supported species. The Fisher statistical test and Yekutieli (FDR under dependency) multitest adjustment methods were applied, and signi cant GO terms with p<0.05 were retrieved. Other parameters were set as the default. All corresponding transcripts were also used in searches against the KEGG database signi cant threshold e-value ≤ 10 −5 . The enriched pathways were downloaded from the KEGG database, and KOBAS (Mao et al., 2005) software was used to test the statistical enrichment of differentially expressed genes in KEGG pathways. Other parameters were set as the default. The functions of lncRNAs were predicted by GO/KEGG enrichment analysis of their predicted target genes.

Interaction network construction
Heat stress-responsive lncRNAs and their target protein-coding genes were included to construct the networks. LncRNAs that were induced or repressed under heat stress in each variety and shared in both were identi ed, and the relationship between them and their target genes was determined. The coexpression network was visualized by Cytoscape software (Shannon et al. 2003).
Colocalization of rice heat stress-related QTLs and heat-responsive lncRNAs QTLs for seedling growth under heat stress were collected from the study of Kilasi et al. (2018). The relevant QTL intervals and their chromosome physical positions were retrieved from the study and compared with the positions of heat-responsive lncRNAs. Those QTLs and lncRNAs with overlapped positions were selected and the QTLs, overlapped lncRNAs and their targets were then mapped on the chromosome with MapChart 2.3 (Voorrips 2002).

Sequence motif search and comparison analysis
For the differentially expressed lncRNAs during high temperature, we packed their body sequences to search the enriched motifs by MEME (Bailey et al. 2009), with the parameters set as default except for motif widths ranging from 6 to 12 nucleotides (Di et al. 2014). To further evaluate the motif signi cance enriched in HRL body sequences, we compared the motif signi cance between HRLs and the background sequences. Three background datasets included upstream 1.5 kb nucleotides of HRLs as predicted HRL promoters and the body sequences of constitutively expressed lncRNAs during heat stress, namely, non-HRL body sequences, as well as 1.5 kb upstream promoter sequences of non-HRL. For motifs identi ed by MEME from HRL body regions, MAST software was applied to search their enrichment in background sequences with an RNA format and p-value less than 0.0001. For each motif in a set of sequences, we calculated the ratio of motif sites to all sequence sites, which was de ned as the length of all sequences divided by the motif width. Then, we compared the statistical signi cance of the motif ratios between different datasets using Fisher's exact test. The -log10(p-value) results from Fisher's test were used to produce a heatmap using an online tool at https://www.omicshare.com/tools/Home/Soft/heatmap.

Results
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 pro ling. For each condition, the seedlings of three biological replications for each of two varieties were used to construct a total of 18 strand-speci c 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-rststrand. 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 Cu inks 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.

Identi cation 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 identi ed. 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-delity lncRNAs were obtained ( Figure 2A).
According to the location of these lncRNAs relative to protein-coding genes, the 2,743 high-delity lncRNAs were classi ed 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 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 Identi cation 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 identi ed 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, cisregulation 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 coe cients (|r| > 0.5 and P < 0.05). Overall, 1278 lncRNA-mRNA co-expression pairs were identi ed, 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 modi cation 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 speci c 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 microRNAgene 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 signi cant 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 identi ed 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 classi ed 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 (speci cally, carbohydrate metabolism, amino-acid metabolism, energy metabolism and so on), genetic information-processing translation, folding, classi cation 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. Thirtyeight 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 speci c 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 coe cients of the steady-state FPKM values between HRLs and the restprotein-encoding genes. Twenty-four lncRNAgene pairs were discovered for heat-responsive lncRNAs shared in both varieties, and 88 and 86 lncRNAgene 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.  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 signi cantly 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 signi cantly 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. Signi cantly 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 signi cantly 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 signi cantly 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 signi cance 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 nger), 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 speci c 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 ndings 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.

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 identi ed 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 heatresponsive lncRNAs identi ed 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 cistarget 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 proteincoding 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 signi cant 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. Elucidating the roles of lncRNAs in regulating rice adaptability to heat stress in this study will provide a novel molecular basis to cope with high temperature. Using steady-state strand-speci c RNA-seq, we identi ed 231 heatresponsive lncRNAs in rice seedlings. An interaction network of HRLs and target-coding genes was constructed to predict the potential roles of HRLs in heat-stress responses. GO enrichment analysis of HRL targets highlighted the general role of lncRNAs in photosynthesis, ATP synthesis, mitochondrialenergy metabolism, transmembrane transporter activity and heat-shock protein Hsp20 at high temperature, independent of variety. We also identi ed the tolerant variety of ZS97B-speci c responsive lncRNAs, which were predicted to potentially interact with coding genes in the nucleoside triphosphate metabolic process, dehydrogenase (NAD+) and Hsp70 protein binding. HRLs in susceptible SYD2 cells were predicted to participate in the protein dimerization activity, apoplast, and RNA polymerase II transcription factor complex pathways. These results re ect the crucial roles of lncRNAs in the main molecular processes triggered by high temperature (Sarkar et  . HSF and Hsp are the most common regulators and factors in response to heat stress. Only three Hsp proteins (Hsp20, Hsp70, and heat-shock protein DnaJ) were predicted to be target genes of heat-responsive lncRNAs. We proposed that HSF and Hsp were involved in the crucial and complex regulatory network under high temperature, and therefore, it is not a simple linear correlation between them and other genes. We expect that more powerful strategies, such as chromatin 3D organization analysis, will provide deep insight into the interaction between  (Ma et al. 2014). We predicted the target-coding genes for all lncRNAs identi ed in this study through the above four methods. Based on the relationship with target genes, we built an interaction network for the heat-responsive lncRNAs in three groups, namely, each accession-speci c HRL and shared HRL in both groups. Except for the above four modes used to predict target genes, the coding genes with absolute Pearson correlation coe cients of FPKM values larger than 0.9 (P < 0.05), but not included in the above four classes, were taken as the trans-target genes and integrated into the network.

Discussion
The cis interaction between HSLs and coding genes represents the most predominant interaction pattern ( Figure 6). Among them, TCONS_00057920 was downregulated during heat stress and modulated 5 neighbouring genes in cis mode. Among these targets, Os03g0190100, encoding prenyltransferase, is involved in the photosynthesis process, and Os03g0190300 encodes a leucine-rich repeat (LRR) family protein, which has been observed to be important in the response to salt stress in Medicago truncatula (de Lorenzo et al. 2009). The important functions of LRR domains in RPP4 in temperature-dependent defence signalling were demonstrated in Arabidopsis (Bao et al. 2014). Other examples, such as TCONS_00100154 and its adjacent downstream gene (Os06g0195800, heat-shock protein DnaJ), were hardly expressed under optimal conditions in both varieties but induced during stress.
LncRNAs have been suggested to serve as competing endogenous RNAs (ceRNAs), which are targeted by speci c miRNAs via target mimics, thereby re ning the interaction between miRNAs and target mRNAs (Franco-Zorrilla et al. 2007;Nejat & Mantri 2018). The ceRNA interaction network among lncRNAs, miRNAs, and mRNAs was constructed in rice under Pi de ciency ). The ceRNA relationships of "HSL-microRNA-gene" were predicted and represented the second largest interacting mode between lncRNA and target genes. For instance, lncRNA TCONS_00043075 ( Figure 6A) is heat responsive in both varieties and as a target mimic of microRNA osa-miR1436. The expression of osa-miR1436 has previously been demonstrated to be induced during heat stress in both tolerant and sensitive genotypes, indicating the general collaborative role of lncRNAs and microRNAs in the heat stress response (Mangrauthia et al. 2017b). Furthermore, 17 protein-coding genes were predicted to form ceRNA relationships with TCONS_00043075 mediated by osa-miR1436. Some of these target genes were involved in important stress-or molecular regulation-relevant processes, such as response to oxidative stress, abiotic stimulus, mRNA capping, protein binding and zinc-ion binding.
HSLs serving as precursors of miRNAs to regulate target-coding genes were predicted. A particularly obvious example is TCONS_00092993, a predicted precursor of osa-miR1850, which was predicted to target 39 protein-coding genes ( Figure 6B). The expression trends of osa-miR1850 and TCONS_00092993 were reversed under heat stress (present study and Mangrauthia et al. 2017a). These target protein-coding genes were enriched in important heat-related processes, such as heat-shock proteins, oxidation-reduction processes, translation, and photosynthesis. Based on these results, we proposed that TCONS_00092993 could play key roles by collaborating with its target genes via osa-miR1850 in plants facing high temperature.
Additionally, transcription factors (TFs) are considered the core regulators in the stress-response process (Wilkins et al. 2016). The centrality of TFs in environmentally in uencing gene networks was highlighted by integrating transcriptome, nucleosome-free chromatin, and cis-regulatory motif occurrence analyses (Wilkins et al. 2016). MicroRNA-mediated TF activity has been reported in rice under drought, heat and salt stresses (Nigam et al. 2015). Among the lncRNA-gene network responsive to heat stress speci cally in SYD2, Gene Os01g0104200 ( Figure 6C) encodes a DNA-binding protein, a subfamily of NAC TFs. Os01g0104200 positively co-expressed with 17 lncRNAs. Since transcription factors could bind DNA and promote gene expression, we speculate that this NAC TF functions as an upstream regulatory factor regulating these lncRNAs, which then regulate other downstream genes with important roles in the heat response in different ways (Results section, Figure 6C).
These results suggest that the interaction network between the HSLs and their target genes highlighted the underlying molecular crosstalk during plant responses to stress. The interaction relationships between HRLs and coding genes predicted in this study could be promising candidates for further validation and future application in plant heat-tolerance improvement.
The combination of genetic mapping and high-throughput sequencing is effective for exploring candidate lncRNAs and coding genes dedicated to heat stress In recent years, diverse and fruitful omics data have grown rapidly. Beyond generating new datasets, another key layer for biological research is to integrate and further explore publicly available data resources. QTL mapping for important agricultural traits is a primary method to investigate the underlying genetic basis regulating phenotype. Some previous studies have reported the genetic variation region affecting the heat-stress response in different rice tissues, such as milky-white grains, seedling Considering the tissue and development stage used in this study, we compared the lncRNAs that were responsive to heat stress with the QTL regions for seedling growth under heat stress (Kilasi et al. 2018). Finally, we found 10 heat-responsive lncRNAs, and some cis-target-coding genes colocalized with QTLs for shoot or root length under heat stress or as a percentage of the control. The cis and ceRNA targets were described as well, and some of them were predicted to have important molecular functions and involvement in heat stress-related biological processes. For example, TCONS_00070624 and its cis-target Os03g0781100 were seated in QTL rlpc3.1, and TCONS_00018800 was adjacent to the cis-target Os10g0368902 located in QTL Slpc10.2. Additionally, although some target genes were not localized within QTLs, such as their relevant lncRNAs, they cooperate via microRNAs to form ceRNA pairs. Heat-responsive lncRNA TCONS_00018499, colocalized with QTL Slpc10.2, can serve as a sponge RNA of osa-miR1439 for 7 coding genes, encoding NADPH HC toxin reductase, protein synthesis and redox. We proposed that the heat-responsive lncRNAs within heat-related QTLs regulated important target-coding genes that are involved in core abiotic stress pathways to cope with a high-temperature environment. The lncRNAs and coding genes identi ed in this section will be important resources for genetic improvement of heat tolerance in rice, such as using them as molecular markers to aid selection.
Speci c motifs are potential target sites for heat-tolerance genetic improvement Previous reports have suggested that RNA molecules with similar functions usually contain conserved sequence motifs (Ray et al. 2013). Since environmentally responsive lncRNAs have stressspeci c expression pro les, stress-related sequence motifs have been enriched in stress-responsive lncRNAs (Di et al. 2014; Zhang et al. 2019). We predicted conserved motifs from heat-responsive lncRNA primary sequences. By comparing the motif enrichment in HRL bodies with background sequences, signi cant 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 were detected. The results indicated the importance of these motifs in speci c lncRNAs functioning in response to heat stress. However, further testing in a greater number of diverse datasets and experimental validation are required to strengthen the ndings here. We anticipate further sequence motif elucidation of abiotic stress-responsive lncRNAs to provide a framework for future biotechnological applications of lncRNAs in crop improvement in abiotic stress tolerance, such as target modi cation relevant to important motifs.

Conclusions
Taken together, to our knowledge, we present the rst comprehensive characterization of heat-responsive lncRNAs (HRLs) in rice. We predicted and constructed the regulated networks among HRLs, coding genes and miRNAs based on the possible mechanisms with which lncRNAs function, such as cis-action, ceRNAs, precursors of miRNA and antisense. The function of HRLs was examined by GO and KEGG enrichment analyses on HRL targets. By integrating HRLs and their targets with publicly available heat-related QTLs, we further strengthened the possible link between these lncRNAs and the heat response in rice. The signi cantly enriched motifs in HRLs will be the best target fragments for further characterization of HRL functions, such as using the CRISPR technique. This study provides new insights into the function of lncRNAs in the heat response and fruitful resources for further investigation on plant heat-tolerance improvement. Further work is needed to experimentally explore in depth the molecular function and mechanisms of the candidate HRLs and target genes in response to high temperature in rice. Additionally, the inadequacies of this study include using the steady-state RNA-sequencing method to measure lncRNA and mRNA expression. Since long noncoding RNAs are unstable, the combination of steady-state and nascent RNA sequencing methods such as global run-on sequencing (GRO-Seq) is necessary in the future to assess RNA transcription and to elucidate the correlation between different RNA species more accurately.  Pipeline for rice lncRNA identi cation and multiomics data integration.    , uniquely in ZS97B upon heat treatment (D) and uniquely in SYD2 upon heat treatment (F). The bubble size represents the number of genes, the colour from green to red indicates the signi cant level of the GO term, the redder colour indicates higher signi cance, and the rich factor indicates the ratio of the "target gene ratio" to the "background gene ratio" enriched in the GO term. A larger rich factor indicates more signi cant enrichment.

Figure 5
Validation of differentially expressed genes and lncRNAs by qRT-PCR. The relative expression levels shown were calculated as the ratio of the investigated genes and lncRNAs to the Actin1 expression level and normalized to 1 in SYD under control conditions. The X-axis indicates the rice materials and growth conditions. The Y axis indicates the relative expression level, and the means ± standard error of n =3 replicates. Actin1 was used as the reference gene. The expression pro les of lncRNAs and relevant protein coding genes are shown at the bottom and top rows, respectively.   The ceRNA pairs between lncRNA TCONS_00004759 and target-coding genes via miRNA osa-miR1438.

Figure 9
The colocalization of heat-responsive lncRNAs with QTLs for seedling growth under heat stress. The blue rectangles indicate the QTL regions, and the red markers indicate lncRNAs within QTLs; the cis-target genes of lncRNAs within QTLs are also marked on the linkage groups and linked with lncRNAs by green arrow lines; the target genes of lncRNAs through ceRNA relationships are not located within QTLs and linked with lncRNAs by red lines of dashes; and the miRNAs of ceRNAs are labelled on the red lines.