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

doi:10.21203/rs.3.rs-740186/v1

Download PDF

Research Article

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

https://doi.org/10.21203/rs.3.rs-740186/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Long noncoding RNAs (lncRNAs) play vital roles in plant responses to environments. However, the systematic identification of lncRNAs in rice under heat stress remains to be fully performed. We performed steady-state strand-specific RNA-seq for rice seedlings under heat stress to predict lncRNAs. In total, 2743 lncRNAs were predicted, and their expression profiles in response to heat treatments were provided. We identified 231 differentially expressed lncRNAs (DELs) in their steady-state level under heat stress, including 31 DELs common to both varieties and 103 and 97 specific to ZS97B and SYD2, respectively, all defined 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 significant motifs enriched in 231 HRL sequences, and the differential significance of enriched motifs between HRLs and background sequences was tested. Our findings provide new insights into the characterization of lncRNAs in terms of the heat response and resources for further investigations of plant heat tolerance improvement.

Medical Genetics

Long noncoding RNAs

rice

heat stress

network

coding gene targets

Identified 2743 rice lncRNAs

LncRNAs in response to heat stress

Function prediction of HRLs

Network among HRLs, genes and miRNAs

Colocalization of HRLs with QTLs

Significant motifs in HRL sequences

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 influences 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 fluidity, reducing the photosynthesis rate, DNA breakage, protein denaturation and so on (Wilkins et al. 2016; Ohama et al. 2017). To some extent, plants can cope with and adapt to environmental changes by coordinating the expression networks of different sets of genes (Nagano et al. 2012; Sarkar et al. 2014). Thus, it is vital to elucidate the molecular mechanisms in plants in response and adaptation to heat stress to promote thermo tolerance improvement (Ohama et al. 2017).

Rice, as one of most important staple foods and a model monocot species, is seriously threatened by rising temperatures from global warming (Peng et al. 2004). The average air temperature is estimated to rise by 1.0-3.7 °C as a result of global climate change by the end of the 21st century (IPCC 2013). Increasing temperatures will dramatically reduce rice production by up to 25% (Wassmann et al. 2009; van Oort and Zwart 2018). Hence, improving rice tolerance to heat stress in terms of sustaining yield stability is of vital importance. Although most relevant studies on rice heat tolerance have mainly focused on the reproductive stage, the impact of heat stress during the vegetative stage also matters by reducing tiller and panicle number mainly through inhibiting photosynthesis, enhancing respiration and oxidative damage and shortening tillering time and growth season (Xu et al. 2020).

In recent years, many important genes and pathways have been identified for their involvement in the response to heat stress in different rice tissues, such as rice seedlings (Mittal et al. 2012), grain at the early milky stage (Liao et al. 2015), spikelets at the meiosis stage (Yang et al. 2015), rice floral organs (Li et al. 2015b) and reproductive tissues during anthesis (Gonzalez-Schain et al. 2016). Co-expression networks suggest that various molecular elements, including heat-shock proteins (Hsps), heat-shock factors (Hsfs), and other transcription factor (TF) families and stress-responsive genes, work together to contribute to rice adaptation to heat stress (Sarkar et al. 2014; (Hu et al. 2009; Mittal et al. 2012; Gonzalez-Schain et al. 2016; Wilkins et al. 2016). Quantitative proteomic analysis has indicated that high-night temperature could disrupt the redox equilibrium of rice cells, triggering calcium-dependent protein kinase and the COP9 signalosome, thereby regulating downstream genes involved in the high-night temperature response (Zhang et al. 2017). Recent data have also identified some genomic regions associated with seedling growth and spikelet fertility in response to high temperature in rice (Kilasi et al. 2018; Nubankoh et al. 2020). GWAS has been conducted to detect quantitative trait loci (QTLs) for rice spikelet sterility at 38 °C, covering genes sensing abiotic stresses and regulating plant responses, such as cell wall-associated kinases, heat-shock proteins and cell division-relevant genes (Lafarge et al. 2017). Combined GWAS and T-DNA insertion knockout line analysis have identified developmental stage-specific QTLs associated with fertility reduction upon heat stress (Bac-Molenaar et al. 2015). Although these efforts have been attempted to elucidate the underlying mechanism of rice responses to heat stress, most key molecules are protein-coding genes. Heat stress-responsive, long noncoding RNAs (lncRNAs) have not been well characterized in rice.

LncRNAs are RNA molecules that are usually untranslated into proteins and longer than 200 bp (Huarte 2015). Emerging evidence strongly suggests the essential roles of lncRNAs in plant development and responses to abiotic stresses at the epigenetic, transcriptional and posttranscriptional levels (Nejat & Mantri 2018). The functions of several lncRNAs have been illustrated in plants involving the regulation of plant flowering, stress response and reproductive growth (Franco-Zorrilla et al. 2007; Wang et al. 2014b; Kim & Sung 2017). 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-specific, 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 identified 1392 spatial- and developmental-specific 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 profiles 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). Identification 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 profiled 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 identified and characterized lncRNAs in two rice varieties with contrasting heat tolerance with the strand-specific RNA-seq (ssRNA-seq) technique. Through analysing the expression profile 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 significant lncRNAs and their key sequence motifs to promote the breeding of rice tolerant to heat stress.

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) filled 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

Total RNA was extracted from the harvested leaves of each cultivar under each treatment by using an E.Z.N.A. plant RNA kit (OMEGA Bio-Tek, USA) and following the manufacturer’s protocol. RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyser 2100 System (Agilent Technologies, CA, USA). A total amount of 1.5 μg RNA per sample was used as input material for rRNA removal using the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, USA). Sequencing libraries were generated using the NEBNextR UltraTM Directional RNA Library Prep Kit for Illumina R (NEB, USA) following the manufacturer’s recommendations. Clustering of the index-coded samples was performed on an acBot Cluster Generation System using TruSeq PE Cluster Kitv3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and paired-end reads were generated.

Read mapping and transcriptome assembly

FastQC analysis was performed to estimate the quality of raw reads. After trimming adapter sequences and filtering low-quality reads with >5% ambiguous bases, the clean reads of each replicate from each cultivar were mapped to the reference genome of Oryza sativa IRGSP-1.0 using the TopHat 2.0.12 program with the default parameter set (Pertea et al. 2016). The transcriptome for each of the 18 libraries was constructed individually using cufflink 2.2.1 with the default parameters. Transcriptome assemblies generated from the above steps were subsequently merged with Cuffmerge to generate comprehensive transcripts for subsequent analysis.

lncRNA identification and characterization

The assembled transcripts were annotated using (i) the Cuffcompare program to filter out the known protein-coding transcripts, rRNA and tRNA with the comparison code‘=’; (ii) transcripts tagged with codes ‘e’, ‘p’ and ‘s’ by Cuffcompare were filtered out, and the remaining transcripts were used to screen for putative lncRNAs; (iii) transcripts with lengths greater than 200 nt and FPKM scores greater than 0.5 in at least one sample were selected as lncRNA candidates; and (iv) further screened through the protein-coding score test using the Coding Potential Calculator (CPC2), Coding-Non-Coding Index (CNCI), Protein family scan software and Protein family database (Pfam) and Plant Long Non-Coding RNA Prediction by Random Forest (PlncPRO) to distinguish the protein-coding genes from the noncoding genes (Sun et al. 2013; Kang et al. 2017; Singh et al. 2017; El-Gebali et al. 2019). Briefly, 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 amino-acid sequences of the intersecting transcripts filtered by CPC2, CNCI and PlncPRO were compared with the pfamA database using PfamScan, and transcripts with E-values smaller than 0.0001 were filtered out. The rest of the transcripts were considered candidate lncRNAs. Based on the location of lncRNAs relative to protein-coding RNAs, lncRNAs were classified 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 identified 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 defined as known lncRNAs.

Gene expression qualification and differential expression analysis

Cufflink 2.2.1 was used to calculate FPKMs of both lncRNAs and coding genes in each sample with default parameters. Gene differential expression analysis was performed using Cuffdiff in Cufflink 2.2.1. Based on the results from Cuffdiff, genes with P-value<0.01, log2FC>1, FDR<0.01 and significant=‘yes’ were assigned as differentially expressed. 18 samples were used for PCA based on gene expression levels with FPKM value bigger than 0.5. Principle component analysis (PCA) and heatmaps generation were performed with https://www.omicshare.com/tools.

qRT-PCR confirmation of DEGs and DELs

Total RNA was extracted from the leaves of each cultivar under each condition using an E.Z.N.A. plant RNA kit (OMEGA Bio-Tek, USA), and three biological replicates were used for qRT-PCR analysis. First-strand cDNA was synthesized by using Superscript III reverse transcriptase (Invitrogen). qRT-PCR was performed on an ABI Prism 7500 real-time PCR system (Applied Biosystems, USA). For qRT-PCR analysis, each lncRNA- and mRNA-specific primer was designed by Primer3 (Table S1). The expression level of the rice Actin1 gene (accession no. X16280) was used as the internal control with primers forward 5′-TGGCATCTCTCAGCACATTCC-3′ and reverse 5′-TGCACAATGGATGGGTCAGA-3′. q-PCR amplifications were performed in an optical 384-well plate. Each reaction was performed in a volume of 25 µL containing 12.5 μL of 2× SYBR green master reagent (Applied Biosystems, COUNTRY NAME), 5.0 μL diluted transcription product, 0.2 μL of each gene-specific primer and 7.1 µL ddH₂O. The following thermal cycle was used: 95 °C for 3 min and then 45 cycles of 95 °C for 30 s, 60 °C for 30 s and 72 °C for 1 min. Dissociation curve analysis was performed using the following thermal profile: 95 °C for 15 s, 60 °C for 20 s and 95 °C for 15 min. The expression levels relative to actin gene were normalized as 2^-△△CT(Livak & Schmittgen 2001).

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 coefficient. 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 identified 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 identified 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 (Xu et al. 2016).

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 classified 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 significant 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 significant threshold e-value ≤ 10⁻⁵. 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 identified, and the relationship between them and their target genes was determined. The co-expression 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 significance enriched in HRL body sequences, we compared the motif significance 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 identified 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 defined as the length of all sequences divided by the motif width. Then, we compared the statistical significance 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.

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.

Extensive involvement of lncRNAs in the response to heat stress in plants revealed by high-throughput sequencing

LncRNAs have been implicated in diverse biological processes, including plant growth, development and environmental stimuli. For instance, lncRNA COOLAIR has been demonstrated to be involved in the epigenetic regulation of FLC during vernalization in response to winter cold (Swiezewski et al. 2009; Kim et al. 2017). LncRNAs were shown to play a possible role in ambient temperature-mediated plant development based upon evidence from COLDWRAP and ASL (Kim & Sung 2017). With the rapid expansion of high-throughput sequencing, a large number of lncRNAs have been identified globally in plants during environmental stresses (Di et al. 2014; Wang et al. 2014a; Shumayla et al. 2017b; Deng et al. 2018; Severing et al. 2018; Wang et al. 2019). As one of the most widespread cereals, rice growth and productivity are severely threatened by global warming. Many studies have been conducted to illustrate the strategies of rice response to stress, but most have focused on coding genes or microRNAs (Kansal et al. 2015; Kilasi et al. 2018; Shi et al. 2018; Nubankoh et al. 2020; Tabassum et al. 2020). 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-specific RNA-seq, we identified 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 identified the tolerant variety of ZS97B-specific 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 reflect the crucial roles of lncRNAs in the main molecular processes triggered by high temperature (Sarkar et al. 2014; Bac-Molenaar et al. 2015; Gonzalez-Schain et al. 2016; Ohama et al. 2017; Zafar et al. 2020). 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 lncRNAs and HSF- or Hsp-coding genes under abiotic stresses (Bohmdorfer & Wierzbicki 2015; Mele & Rinn 2016).

Heat-responsive lncRNAs play regulatory functions through multiple interaction modes with coding genes

Previous studies have suggested that lncRNAs can perform their function through diverse mechanisms, such as regulating the expression of their neighbouring genes by cis action (SI et al. 2018), acting as sponges for miRNA binding to become ceRNAs for the target gene (Cesana et al. 2011), modulating the sense-coding gene as an antisense transcript (Liu et al. 2018b) or as precursors of microRNAs (Ma et al. 2014). We predicted the target-coding genes for all lncRNAs identified 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-specific 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 coefficients 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. 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 specific miRNAs via target mimics, thereby refining 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 deficiency (Xu et al. 2016). 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 influencing 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 specifically 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 growth, and spikelets (Miyahara et al. 2017; Kilasi et al. 2018; Nubankoh et al. 2020). 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 identified 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.

Specific 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 stress-specific expression profiles, 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, 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 were detected. The results indicated the importance of these motifs in specific 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 findings 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 modification relevant to important motifs.

Taken together, to our knowledge, we present the first 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 significantly 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.

Author contributions: B. X. designed the research; Z.Z., H.Z. and B.N. analyzed data; H.Z., B.X. and B.N. conduct the experiments; Z.Z and B.X. wrote the manuscript.

Conflict of interest

No conflict of interest declared.

Acknowledgments

This work was supported by the self-determined research fund of Central China Normal University (CCNU18QN027) and the National Special Key Project of China on Transgenic Research (2016ZX 08001-003).

Bac-Molenaar J.A., Fradin E.F., Becker F.F., Rienstra J.A., van der Schoot J., Vreugdenhil D. & Keurentjes J.J. (2015) Genome-Wide Association Mapping of Fertility Reduction upon Heat Stress Reveals Developmental Stage-Specific QTLs in Arabidopsis thaliana. Plant Cell 27, 1857-74.
Bailey T.L., Boden M., Buske F.A., Frith M., Grant C.E., Clementi L., Ren J., Li W.W. & Noble W.S. (2009) MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res 37, W202-8.
Bao F., Huang X., Zhu C., Zhang X., Li X. & Yang S. (2014) Arabidopsis HSP90 protein modulates RPP4-mediated temperature-dependent cell death and defense responses. New Phytol 202, 1320-34.
Bardou F., Ariel F., Simpson C., Romero-Barrios N., Laporte P., Balzergue S., Brown J. & Crespi M. (2014) - Long noncoding RNA modulates alternative splicing regulators in Arabidopsis. Dev Cell 30, 166-76.
Bohmdorfer G. & Wierzbicki A.T. (2015) Control of Chromatin Structure by Long Noncoding RNA. Trends Cell Biol 25, 623-32.
Cagirici H., Alptekin B. & Budak H. (2017) RNA Sequencing and Co-expressed Long Non-coding RNA in Modern and Wild Wheats. Sci Rep 7, 017-11170.
Cesana M., Cacchiarelli D., Legnini I., Santini T., Sthandier O., Chinappi M., Tramontano A. & Bozzoni I. (2011) A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 147, 358-69.
Chekanova J.A. (2015) Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol 27, 207-16.
Csorba T., Questa J., Sun Q. & Dean C. (2014) Antisense COOLAIR mediates the coordinated switching of chromatin states at FLC during vernalization. Proc Natl Acad Sci U S A 111, 16160-5.
Dai X., Zhuang Z. & Zhao P.X. (2018) psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Res 46, W49-W54.
de Lorenzo L., Merchan F., Laporte P., Thompson R., Clarke J., Sousa C. & Crespi M. (2009) A novel plant leucine-rich repeat receptor kinase regulates the response of Medicago truncatula roots to salt stress. Plant Cell 21, 668-80.
Deng F., Zhang X., Wang W., Yuan R. & Shen F. (2018) Identification of Gossypium hirsutum long non-coding RNAs (lncRNAs) under salt stress. BMC Plant Biol 18, 23.
Di C., Yuan J., Wu Y., Li J., Lin H., Hu L., Zhang T., Qi Y., Gerstein M.B., Guo Y. & Lu Z.J. (2014) Characterization of stress-responsive lncRNAs in Arabidopsis thaliana by integrating expression, epigenetic and structural features. Plant J 80, 848-61.
Ding J., Lu Q., Ouyang Y., Mao H., Zhang P., Yao J., Xu C., Li X., Xiao J. & Zhang Q. (2012) A long noncoding RNA regulates photoperiod-sensitive male sterility, an essential component of hybrid rice. Proc Natl Acad Sci U S A 109, 2654-9.
El-Gebali S., Mistry J., Bateman A., Eddy S.R., Luciani A., Potter S.C., Qureshi M., Richardson L.J., Salazar G.A., Smart A., Sonnhammer E.L.L., Hirsh L., Paladin L., Piovesan D., Tosatto S.C.E. & Finn R.D. (2019) The Pfam protein families database in 2019. Nucleic Acids Res 47, D427-D32.
Franco-Zorrilla J.M., Valli A., Todesco M., Mateos I., Puga M.I., Rubio-Somoza I., Leyva A., Weigel D., Garcia J.A. & Paz-Ares J. (2007) Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet 39, 1033-7.
Gai Y., Yuan S., Zhao Y., Zhao H., Zhang H. & Ji X. (2018) A Novel LncRNA, MuLnc1, Associated With Environmental Stress in Mulberry (Morus multicaulis). Front Plant Sci 9.
Gonzalez-Schain N., Dreni L., Lawas L., Galbiati M., Colombo L., Heuer S., Jagadish K. & Kater M. (2016) Genome-Wide Transcriptome Analysis During Anthesis Reveals New Insights into the Molecular Basis of Heat Stress Responses in Tolerant and Sensitive Rice Varieties. Plant Cell Physiol 57, 57-68.
Hsu K.H., Liu C.C., Wu S.J., Kuo Y.Y., Lu C.A., Wu C.R., Lian P.J., Hong C.Y., Ke Y.T., Huang J.H. & Yeh C.H. (2014) Expression of a gene encoding a rice RING zinc-finger protein, OsRZFP34, enhances stomata opening. Plant Mol Biol 86, 125-37.
Hu W., Hu G. & Han B. (2009) Genome-wide survey and expression profiling of heat shock proteins and heat shock factors revealed overlapped and stress specific response under abiotic stresses in rice. Plant Sci 176, 583-90.
Huarte M. (2015) The emerging role of lncRNAs in cancer. Nat Med 21, 1253-61.
IPCC (2013). Working Group I Contribution to the IPCC Fifth Assessment Report on Climate Change: The Physical Science Basis, Summary for Policymakers.
Jin J., Tian F., Yang D.C., Meng Y.Q., Kong L., Luo J. & Gao G. (2017) PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res 45, D1040-D5.
Kang Y.J., Yang D.C., Kong L., Hou M., Meng Y.Q., Wei L. & Gao G. (2017) CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res 45, W12-W6.
Kansal S., Mutum R., Balyan S.C., Arora M.K., Singh A.K., Mathur S. & Raghuvanshi S. (2015) Unique miRNome during anthesis in drought-tolerant indica rice var. Nagina 22. Planta 241, 1543-59.
Kilasi N.L., Singh J., Vallejos C.E., Ye C., Jagadish S.V.K., Kusolwa P. & Rathinasabapathi B. (2018) Heat Stress Tolerance in Rice (Oryza sativa L.): Identification of Quantitative Trait Loci and Candidate Genes for Seedling Growth Under Heat Stress. Front Plant Sci 9, 1578.
Kim D.H. & Sung S. (2017) Vernalization-Triggered Intragenic Chromatin Loop Formation by Long Noncoding RNAs. Dev Cell 40, 302-12 e4.
Kim D.H., Xi Y. & Sung S. (2017) Modular function of long noncoding RNA, COLDAIR, in the vernalization response. PLoS Genet 13, e1006939.
Kobayashi K., Kawabata M., Hisano K., Kazama T., Matsuoka K., Sugita M. & Nakamura T. (2012) Identification and characterization of the RNA binding surface of the pentatricopeptide repeat protein. Nucleic Acids Res 40, 2712-23.
Lafarge T., Bueno C., Frouin J., Jacquin L., Courtois B. & Ahmadi N. (2017) Genome-wide association analysis for heat tolerance at flowering detected a large set of genes involved in adaptation to thermal and other stresses. PLoS One 12, e0171254.
Lancaster L.T. & Humphreys A.M. (2020) Global variation in the thermal tolerances of plants. Proc Natl Acad Sci U S A.
Li J., Ma W., Zeng P., Wang J., Geng B., Yang J. & Cui Q. (2015a) LncTar: a tool for predicting the RNA targets of long noncoding RNAs. Brief Bioinform 16, 806-12.
Li J.F., Chung H.S., Niu Y., Bush J., McCormack M. & Sheen J. (2013) Comprehensive protein-based artificial microRNA screens for effective gene silencing in plants. Plant Cell 25, 1507-22.
Li L., Eichten S.R., Shimizu R., Petsch K., Yeh C.T., Wu W., Chettoor A.M., Givan S.A., Cole R.A., Fowler J.E., Evans M.M., Scanlon M.J., Yu J., Schnable P.S., Timmermans M.C., Springer N.M. & Muehlbauer G.J. (2014) Genome-wide discovery and characterization of maize long non-coding RNAs. Genome Biol 15, R40.
Li X., Lawas L.M., Malo R., Glaubitz U., Erban A., Mauleon R., Heuer S., Zuther E., Kopka J., Hincha D.K. & Jagadish K.S. (2015b) Metabolic and transcriptomic signatures of rice floral organs reveal sugar starvation as a factor in reproductive failure under heat and drought stress. Plant Cell Environ 38, 2171-92.
Liao J., Zhou H., Peng Q., Zhong P., Zhang H., He C. & Huang Y. (2015) Transcriptome changes in rice (Oryza sativa L.) in response to high night temperature stress at the early milky stage. BMC Genomics 16, 015-1222.
Liu J., Jung C., Xu J., Wang H., Deng S., Bernad L., Arenas-Huertero C. & Chua N.-H. (2012) Genome-Wide Analysis Uncovers Regulation of Long Intergenic Noncoding RNAs in <em>Arabidopsis</em>. The Plant Cell 24, 4333-45.
Liu W., C C., Y L., X X., Z L. & http://orcid.org/---621X I.-O. (2018a) - Genome-wide identification and characterization of mRNAs and lncRNAs involved in cold stress in the wild banana (Musa itinerans). PLoS One 13.
Liu X., Li D., Zhang D., Yin D., Zhao Y., Ji C., Zhao X., Li X., He Q., Chen R., Hu S. & Zhu L. (2018b) A novel antisense long noncoding RNA, TWISTED LEAF, maintains leaf blade flattening by regulating its associated sense R2R3-MYB gene in rice. New Phytol 218, 774-88.
Livak K.J. & Schmittgen T.D. (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego, Calif.) 25, 402-8.
Luo S., Lu J.Y., Liu L., Yin Y., Chen C., Han X., Wu B., Xu R., Liu W., Yan P., Shao W., Lu Z., Li H., Na J., Tang F., Wang J., Zhang Y.E. & Shen X. (2016) Divergent lncRNAs Regulate Gene Expression and Lineage Differentiation in Pluripotent Cells. Cell Stem Cell 18, 637-52.
Ma X., Shao C., Jin Y., Wang H. & Meng Y. (2014) Long non-coding RNAs: a novel endogenous source for the generation of Dicer-like 1-dependent small RNAs in Arabidopsis thaliana. RNA Biol 11, 373-90.
Mangrauthia S., Bhogireddy S., Agarwal S., Prasanth V., Voleti S., Neelamraju S. & Subrahmanyam D. (2017a) Genome-wide changes in microRNA expression during short and prolonged heat stress and recovery in contrasting rice cultivars. J Exp Bot 68, 2399-412.
Mangrauthia S.K., Bhogireddy S., Agarwal S., Prasanth V.V., Voleti S.R., Neelamraju S. & Subrahmanyam D. (2017b) Genome-wide changes in microRNA expression during short and prolonged heat stress and recovery in contrasting rice cultivars. J Exp Bot 68, 2399-412.
Mele M. & Rinn J.L. (2016) "Cat's Cradling" the 3D Genome by the Act of LncRNA Transcription. Mol Cell 62, 657-64.
Mittal D., Madhyastha D.A. & Grover A. (2012) Genome-wide transcriptional profiles during temperature and oxidative stress reveal coordinated expression patterns and overlapping regulons in rice. PLoS One 7, e40899.
Miyahara K., Wada T., Sonoda J.Y., Tsukaguchi T., Miyazaki M., Tsubone M., Yamaguchi O., Ishibashi M., Iwasawa N., Umemoto T. & Kondo M. (2017) Detection and validation of QTLs for milky-white grains caused by high temperature during the ripening period in Japonica rice. Breed Sci 67, 333-9.
Mueller S.P., Unger M., Guender L., Fekete A. & Mueller M.J. (2017) Phospholipid:Diacylglycerol Acyltransferase-Mediated Triacylglyerol Synthesis Augments Basal Thermotolerance. Plant Physiol 175, 486-97.
Nagano A.J., Sato Y., Mihara M., Antonio B.A., Motoyama R., Itoh H., Nagamura Y. & Izawa T. (2012) Deciphering and prediction of transcriptome dynamics under fluctuating field conditions. Cell 151, 1358-69.
Nejat N. & Mantri N. (2018) Emerging roles of long non-coding RNAs in plant response to biotic and abiotic stresses. Crit Rev Biotechnol 38, 93-105.
Nigam D., Kumar S., Mishra D., Rai A., Smita S. & Saha A. (2015) - Synergistic regulatory networks mediated by microRNAs and transcription factors under drought, heat and salt stresses in Oryza Sativa spp. Gene 555, 127-39.
Nubankoh P., Wanchana S., Saensuk C., Ruanjaichon V., Cheabu S., Vanavichit A., Toojinda T., Malumpong C. & Arikit S. (2020) QTL-seq reveals genomic regions associated with spikelet fertility in response to a high temperature in rice (Oryza sativa L.). Plant Cell Rep 39, 149-62.
Ohama N., Sato H., Shinozaki K. & Yamaguchi-Shinozaki K. (2017) Transcriptional Regulatory Network of Plant Heat Stress Response. Trends Plant Sci 22, 53-65.
Paytuvi Gallart A., Hermoso Pulido A., Anzar Martinez de Lagran I., Sanseverino W. & Aiese Cigliano R. (2016) GREENC: a Wiki-based database of plant lncRNAs. Nucleic Acids Res 44, D1161-6.
Peng S., Huang J., Sheehy J.E., Laza R.C., Visperas R.M., Zhong X., Centeno G.S., Khush G.S. & Cassman K.G. (2004) Rice yields decline with higher night temperature from global warming. Proc Natl Acad Sci U S A 101, 9971-5.
Pertea M., Kim D., Pertea G.M., Leek J.T. & Salzberg S.L. (2016) Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc 11, 1650-67.
Qi X., Xie S., Liu Y., Yi F. & Yu J. (2013) Genome-wide annotation of genes and noncoding RNAs of foxtail millet in response to simulated drought stress by deep sequencing. Plant Mol Biol 83, 459-73.
Qin T., Zhao H., Cui P., Albesher N. & Xiong L. (2017) A Nucleus-Localized Long Non-Coding RNA Enhances Drought and Salt Stress Tolerance. Plant Physiol 175, 1321-36.
Ray D., Kazan H., Cook K.B., Weirauch M.T., Najafabadi H.S., Li X., Gueroussov S., Albu M., Zheng H., Yang A., Na H., Irimia M., Matzat L.H., Dale R.K., Smith S.A., Yarosh C.A., Kelly S.M., Nabet B., Mecenas D., Li W., Laishram R.S., Qiao M., Lipshitz H.D., Piano F., Corbett A.H., Carstens R.P., Frey B.J., Anderson R.A., Lynch K.W., Penalva L.O., Lei E.P., Fraser A.G., Blencowe B.J., Morris Q.D. & Hughes T.R. (2013) A compendium of RNA-binding motifs for decoding gene regulation. Nature 499, 172-7.
Sarkar N.K., Kim Y.K. & Grover A. (2014) Coexpression network analysis associated with call of rice seedlings for encountering heat stress. Plant Mol Biol 84, 125-43.
Severing E., Faino L., Jamge S., Busscher M., Kuijer-Zhang Y., Bellinazzo F., Busscher-Lange J., Fernandez V., Angenent G.C., Immink R.G.H. & Pajoro A. (2018) Arabidopsis thaliana ambient temperature responsive lncRNAs. BMC Plant Biol 18, 145.
Shannon P., Markiel A., Ozier O., Baliga N.S., Wang J.T., Ramage D., Amin N., Schwikowski B. & Ideker T. (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13, 2498-504.
Shi W., Li X., Schmidt R.C., Struik P.C., Yin X. & Jagadish S.V.K. (2018) Pollen germination and in vivo fertilization in response to high-temperature during flowering in hybrid and inbred rice. Plant Cell Environ 41, 1287-97.
Shuai P., Liang D., Tang S., Zhang Z., Ye C., Su Y., Xia X. & Yin W. (2014) - Genome-wide identification and functional prediction of novel and drought-responsive lincRNAs in Populus trichocarpa. J Exp Bot 65, 4975-83.
Shumayla, S S., M T., S T., K S. & SK U. (2017a) - Survey of High Throughput RNA-Seq Data Reveals Potential Roles for lncRNAs during Development and Stress Response in Bread Wheat. Front Plant Sci 8.
Shumayla, Sharma S., Taneja M., Tyagi S., Singh K. & Upadhyay S.K. (2017b) Survey of High Throughput RNA-Seq Data Reveals Potential Roles for lncRNAs during Development and Stress Response in Bread Wheat. Front Plant Sci 8, 1019.
SI A., L Z., BR G. & Id O. (2018) - Vex-seq: high-throughput identification of the impact of genetic variation on pre-mRNA splicing efficiency. Genome Biol 19, 018-1437.
Singh U., Khemka N., Rajkumar M., Garg R. & Jain M. (2017) PLncPRO for prediction of long non-coding RNAs (lncRNAs) in plants and its application for discovery of abiotic stress-responsive lncRNAs in rice and chickpea. Nucleic Acids Res 45.
Song X., Liu G., Huang Z., Duan W., Tan H., Li Y. & Hou X. (2016) Temperature expression patterns of genes and their coexpression with LncRNAs revealed by RNA-Seq in non-heading Chinese cabbage. BMC Genomics 17, 297.
Stirpe F., Barbieri L., Gorini P., Valbonesi P., Bolognesi A. & Polito L. (1996) Activities associated with the presence of ribosome-inactivating proteins increase in senescent and stressed leaves. FEBS Lett 382, 309-12.
Sun L., Luo H., Bu D., Zhao G., Yu K., Zhang C., Liu Y., Chen R. & Zhao Y. (2013) Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res 41, 27.
Swiezewski S., Liu F., Magusin A. & Dean C. (2009) Cold-induced silencing by long antisense transcripts of an Arabidopsis Polycomb target. Nature 462, 799-802.
Szczesniak M.W., Rosikiewicz W. & Makalowska I. (2016) CANTATAdb: A Collection of Plant Long Non-Coding RNAs. Plant Cell Physiol 57, e8.
Tabassum R., Dosaka T., Ichida H., Morita R., Ding Y., Abe T. & Katsube-Tanaka T. (2020) FLOURY ENDOSPERM11-2 encodes plastid HSP70-2 involved with the temperature-dependent chalkiness of rice (Oryza sativa L.) grains. Plant J.
Thomson D.W. & Dinger M.E. (2016) Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet 17, 272-83.
van Oort, P. A. J., and Zwart, S. J. (2018). Impacts of climate change on rice production in Africa and causes of simulated yield changes. Glob. Chang. Biol. 24, 1029–1045.
Wang A., J H., C G., G C., B W., C L., L S., Y D. & G Z. (2019) - Genome-wide analysis of long non-coding RNAs unveils the regulatory roles in the heat tolerance of Chinese cabbage (Brassica rapa ssp.chinensis). Sci Rep 9, 019-41428.
Wang H., Liu J., Jang I.-C., Kean M.J., Xu J. & Chua N.-H. (2014a) Genome-wide identification of long noncoding natural antisense transcripts and their responses to light in Arabidopsis. Genome Res 24, 444-53.
Wang M., Yuan D., Tu L., Gao W., He Y., Hu H., Wang P., Liu N., Lindsey K. & Zhang X. (2015) Long noncoding RNAs and their proposed functions in fibre development of cotton (Gossypium spp.). New Phytol 207, 1181-97.
Wang M., Zhao W., Gao L. & Zhao L. (2018) Genome-wide profiling of long non-coding RNAs from tomato and a comparison with mRNAs associated with the regulation of fruit ripening. BMC Plant Biol 18, 75.
Wang Y., Fan X., Lin F., He G., Terzaghi W., Zhu D. & Deng X.W. (2014b) Arabidopsis noncoding RNA mediates control of photomorphogenesis by red light. Proc Natl Acad Sci U S A 111, 10359-64.
Wassmann, R., Jagadish, S. V. K., Heuer, S., Ismail, A., Redona, E., Serraj, R., et al. (2009). Climate change affecting rice production: the physiological and agronomic basis for possible adaptation strategies. Adv. Agron. 101, 59–122.
Wilkins O., Hafemeister C., Plessis A., Holloway-Phillips M.M., Pham G.M., Nicotra A.B., Gregorio G.B., Jagadish S.V., Septiningsih E.M., Bonneau R. & Purugganan M. (2016) EGRINs (Environmental Gene Regulatory Influence Networks) in Rice That Function in the Response to Water Deficit, High Temperature, and Agricultural Environments. Plant Cell 28, 2365-84.
Wilusz J.E., Sunwoo H. & Spector D.L. (2009) Long noncoding RNAs: functional surprises from the RNA world. Genes Dev 23, 1494-504.
Xu, J., A. Henry, et al. (2020). "Rice yield formation under high day and night temperatures-A prerequisite to ensure future food security." Plant Cell Environ 43(7): 1595-1608
Xu X.W., Zhou X.H., Wang R.R., Peng W.L., An Y. & Chen L.L. (2016) Functional analysis of long intergenic non-coding RNAs in phosphate-starved rice using competing endogenous RNA network. Sci Rep 6, 20715.
Yang J., Chen X., Zhu C., Peng X., He X., Fu J., Ouyang L., Bian J., Hu L., Sun X., Xu J. & He H. (2015) RNA-seq reveals differentially expressed genes of rice (Oryza sativa) spikelet in response to temperature interacting with nitrogen at meiosis stage. BMC Genomics 16, 959.
Yin P., Li Q., Yan C., Liu Y., Liu J., Yu F., Wang Z., Long J., He J., Wang H.W., Wang J., Zhu J.K., Shi Y. & Yan N. (2013) Structural basis for the modular recognition of single-stranded RNA by PPR proteins. Nature 504, 168-71.
Zafar S.A., Hameed A., Ashraf M., Khan A.S., Qamar Z.U., Li X. & Siddique K.H.M. (2020) Agronomic, physiological and molecular characterisation of rice mutants revealed the key role of reactive oxygen species and catalase in high-temperature stress tolerance. Funct Plant Biol 47, 440-53.
Zhang H.Y., Lei G., Zhou H.W., He C., Liao J.L. & Huang Y.J. (2017) Quantitative iTRAQ-based proteomic analysis of rice grains to assess high night temperature stress. Proteomics 17.
Zhang Z., Zheng Y., Ham B.K., Zhang S., Fei Z. & Lucas W.J. (2019) Plant lncRNAs are enriched in and move systemically through the phloem in response to phosphate deficiency. J Integr Plant Biol 61, 492-508.

Figuresinmanuscriptfinalv2supplementsubmitpicture.pptx
Figure S1. Seedlings growth of ZS97B and SYD2 under heat treatment at 0, 2 and 6 days. Figure S2. Principal component analysis of samples in this study. a, b, c indicates three replicates of each variety under a specific condition. Figure S3. Interaction network of ceRNAs(lncRNAs) identified in this study with their targets and relevant miRNAs. Green circle, lncRNA; red circle, miRNAs; blue circle, genes. Figure S4. GO and KEGG enrichment analysis of coding gene targets of lncRNAs identified in this study. (A) GO enrichment analysis of coding gene targets of lncRNAs identified in this study. (B) KEGG enrichment analysis of lncRNA targets. Figure S5 The KEGG oxidative phosphorylation pathway including coding gene targets of DELs under heat treatments. The KEGG Ribosome pathway map can be found online at http://www.kegg.jp/pathway. The red boxes indicate the proteins encoded by the gene targets of DELs under heat treatments. Figure S6. The KEGG Ribosome pathway including coding gene targets of DELs under heat treatments. The KEGG Ribosome pathway map can be found online at http://www.kegg.jp/pathway. The red boxes indicate the proteins encoded by the gene targets of DELs under heat treatments. Figure S7. The KEGG photosynthesis pathway including coding gene targets of DELs under heat treatments. The KEGG Ribosome pathway map can be found online at http://www.kegg.jp/pathway. The red boxes indicate the proteins encoded by the gene targets of DELs under heat treatments. Figure S8. The KEGG vitamin B 6 metabolism pathway including coding gene targets of DELs under heat treatments. The KEGG Ribosome pathway map can be found online at http://www.kegg.jp/pathway. The red boxes indicate the proteins encoded by the gene targets of DELs under heat treatments. Figure S9. Interaction network of lncRNAs differentially expressed in both varieties upon heat treatment and their target coding genes and relevant miRNAs. The relationship of lncRNAs and mRNAs shown in different colors or shapes of lines as indicated in the legends. Red ellipse, lncRNA; green circle, mRNAs; yellow Octagon, miRNAs. Figure S10. Interaction network of lncRNAs differentially expressed in ZS97B upon heat treatment and their target coding genes and relevant miRNAs. The relationship of lncRNAs and mRNAs shown in different colors or shapes of lines as indicated in the legends in Figure S7. Red ellipse, lncRNA; green circle, mRNAs; yellow Octagon, miRNAs. Figure S11. Interaction network of lncRNAs differentially expressed in SYD2 upon heat treatment and their target coding genes and relevant miRNAs. The relationship of lncRNAs and mRNAs shown in different colors or shapes of lines as indicated in the legends in Figure S7. Red ellipse, lncRNA; green circle, mRNAs; yellow Octagon, miRNAs.
ManuscriptsSupplementalTablesv2.xlsx
Table S1 Primers used for qRT-PCR. Table S2 RNA-seq data for 18 samples. Table S3 Statistics for Alignment of the Clean Reads on Reference Genome. Table S4 cis target genes of lncRNAs. Table S5 Antisense lncRNA and their target mRNA. Table S6 ceRNAs pair of lncRNA and mRNA. Table S7 lncRNA(as miRNA precursor)-miRNA-target-gene. Table S8 DELs shared in both varities interactions with target coding genes. Table S9 DELs in ZS97B interactions with target coding genes. Table S10 DELs in SYD interactions with target coding genes. Table S11 KEGG enrichment analysis of heat responsive lncRNAs by their target genes.

Download PDF

Reviewers invited by journal
31 Aug, 2021
Reviews received at journal
05 Aug, 2021
Editor assigned by journal
22 Jul, 2021
First submitted to journal
21 Jul, 2021

You are reading this latest preprint version

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

Status:

Version 1

Abstract

Figures

Key Message

Introduction

Methods

Results

Discussion

Conclusions

Declarations

References

Supplementary Files

Status:

Version 1