Transcriptional Landscape of Pathogen-responsive Lncrnas Reveals OS_LNC1812-Induced Disease Resistance by Abscisic Acid Pathway in Rice

Background: Long non-coding RNAs (lncRNAs) produced by the plant genome are essential regulators in diverse biological processes. Despite increasing knowledge of lncRNAs in plant development, little is known about responding to plant disease, especially in rice. Results: In this study, we investigated a comprehensive disease-responding lncRNA proles in rice defence against Rice black-streaked dwarf virus (RBSDV) and Rice stripe virus (RSV) for the rst time. Analysis of RNA sequencing of rice leaves infected with the two pathogens identied total 1,925 lncRNAs, in which 724 lncRNA were derived from alternative splicing (AS) events and the largest number of lncRNAs produced by intron retention (IR) of AS events. We calculated the regulatory relationship between transcription factors, lncRNAs, and mRNAs by constructing gene regulatory network and found that a lncRNA, OS_LNC1812 regulated by Os01t0730700-01 encoding transcription factor WRKY14 may target abscisic acid (ABA) responsive element binding factor ABF (Os01t0867300-01) involved in ABA signaling pathway during pathogens infection. We subsequently found that highly expression of OS_LNC1812, WRKY14, and ABF after ABA treatment. Conclusions: Our results reveal the common lncRNAs proles respond to RBSDV and RSV in rice and provide novel potential modulation of lncRNAs in ABA pathway in the regulation of plant disease resistance.


Background
As one of the most important cereal crops, rice (Oryza sativa L.) is widely consumed by more than half of population in the world. However, multiple diseases severely decrease rice yield and quality, such as rice stripe virus disease (RSVD) and rice black-streaked dwarf disease (RBSDD). RSVD, one of the most serious viral diseases to rice, is caused by rice stripe virus (RSV). In addition to rice, other staple crops are also infected by RSV, such as barley, wheat, and maize. After RSVD, RBSDD caused by rice black-streaked dwarf virus (RBSDV) is also an important disease of rice from many Asian countries, such as China, and Japan [1]. Similar to RSV, RBSDV can also harm other staple crops, such as wheat (wheat dark-green dwarf disease), maize (Maize rough dwarf disease) [2,3], and several species of weeds [4][5][6][7]. The most effective, economical and safe strategy to control RSV and RBSDV diseases is breeding resistant cultivars [8,9].
Investigating resistance mechanisms by identifying resistance genes is an important contributor for breeding resistant cultivars. To data, many resistance genes were identi ed in rice response to RSV or RBSDV infection. Overexpression of OsCIPK30 can enhance the tolerance in rice to RSV by positively regulating pathogenesis-related genes [10]. Ubiquitin like protein 5 participates in the resistance mechanism to RSV infection by mediating the degradation of RSV p3 protein via the 26S proteasome in rice and Nicotiana benthamiana [11]. During RBSDV infection, The up-regulated rice signal transductionassociated genes can improve rice resistance by inducting of the jasmonic acid (JA) pathway and inhibiting the abscisic acid (ABA) pathway [12,13]. Recently reports shown that in addition to resistance genes, non-coding RNAs (ncRNAs) including small ncRNAs and long ncRNAs (lncRNAs) also play an important in rice defense against pathogens [14,15].
Although many works associated with functional characterization of small ncRNAs are reported, the function of lncRNAs in plants is largely unknown. In plant, transcriptome-wide studies shown that lots of lncRNAs are transcribed from genome, which is an important contributor in almost all biological processes [16,17]. Generally, transcript with no protein-coding ability and length > 200bp is de ned as lncRNA [18]. According to the genomic location, lncRNAs are divided into three major categories, long intronic noncoding RNA, long noncoding natural antisense transcripts (lncNAT), and long intergenic noncoding RNA (lincRNA) [19]. There are many reports about lncRNAs in multiple plant species [20], such as Solanum lycopersicum [21], Zea mays [22], Oryza sativa [23], Arachis hypogaea [24,25], Populus [26] and Arabidopsis thaliana [27]. However, most of studies on lncRNAs focus on plant growth, development and abiotic stress responses, such as growth [28,29], fruit ripening [30], salt and drought stresses [31], and nutritional stimulation [32]. In contrast, lncRNAs involved in resistance regulatory mechanisms against pathogens are rarely reported, especially in rice [33,34].
Tianze et al revealed the regulatory network between lncRNA and mRNA in rice response to RBSDV infection by performing RNA-seq [35]. To date, the expression pro le and functional characterization of lncRNAs during RSV and RBSDV infection in rice are still largely unknown. Therefore, to investigate the rice immune defense against RSV and RBSDV infection, we genome-wide identi ed rice lncRNAs and investigated the common regulatory mechanism of rice response to RSV and RBSDV infection for the rst time. Our results provided the valuable information for investigating defense mechanisms and breeding resistant rice.

Results
Global identi cation of lncRNAs in Oryza sativa L To globally identity lncRNAs associated with response to RBSDV and RSV infection in rice, we used RBSDV and RSV to inoculate rice respectively for 30 days and then performed RNA-seq using diseased rice. There were 9 cDNA libraries constructed from rice treated with RBSDV, RSV, and CK (control) (3 biological replicates for each sample). We performed RT-PCR and found that RBSDV and RSV were successfully colonized in corresponding samples (Additional le 1: Fig. S1; Additional le 5: Table S1).
We used Illumina HiSeq 4000 platform to sequence all of libraries and obtained a total of 0.62 billion raw reads. There were 0.56 billion clean reads generated after removing adaptor sequences and low-quality reads. The average percentage of clean reads and GC content were 90.48% and 48.51%, respectively (Additional le 6: Table S2). Then, approximately 78.46%-96.08% of clean reads were mapped to Oryza sativa L. reference genome (IRGSP-1.0) by performing HISAT2 software (Additional le 7: Table S3). Total 46,020 transcripts (38,819 genes) including 24,893 known transcripts and 21,127 novel transcripts were generated. We calculated the Transcripts Per Million (TPM) represented the expression level of all transcripts by using Salmon software. We considered 41,992 transcripts (20,240 mRNAs) with TPM > 1 at least one sample to be expression (Additional le 8: Table S4). We used the novel transcripts with TPM > 1 at least one sample to identify lncRNAs by the ltering pipeline (Fig. 1a). The transcripts with length < 200bp, number of exon = 1, and ORF > 300bp were ltered. Subsequently, CPC software, and Swiss-Prot database were used to evaluate protein-coding potential. Finally, we obtained a total of 1,925 lncRNAs including 1,112 genic lncRNAs, and 813 intergenic lncRNAs ( Fig. 1b; Additional le 9: Table S5). In these lncRNAs, there were 752, 418, 395, and 360 sense-genic, sense-intergenic, antisense-intergenic, and antisense-genic lncRNAs, respectively (Fig. 1c). Sequence lengths between 200 and 2000 bp of lncRNAs and mRNAs were most frequent and their percentages reached 72.49% and 70.75%, respectively (Fig. 1d).

Characterization of rice lncRNAs
To analysis the distribution on rice chromosomes, we performed the R package Circlize() for visualization of the location of lncRNAs and mRNAs on rice chromosomes. The result indicated that lncRNAs were evenly distributed on 12 rice chromosomes, while higher densities of mRNA was located in the chromosome "arms" of most rice chromosomes than pericentromeric regions (Fig. 2a). Expression analysis showed that the expression level of identi ed lncRNAs was lower than mRNAs in rice (Fig. 2b). The different distribution tendencies of exon number were observed and approximately 2,731 (14.43%) mRNAs and 770 (40%) lncRNAs were composed of two exons (Fig. 2c). AS events, a way to generate mature transcripts by excising introns, was globally identi ed by performing AStalavista software. We obtained total 23,898 AS events involved in 5,774 genes including 724 lncRNAs and 4,424 mRNAs ( Fig.   2d; Additional le 10: Table S6). Subsequently, we customized a user-friendly program to identify the ve major AS events, AA, AD, IR, ES, and MX. IR was the most abundant (6,486,27.14%) of all AS events and AA was followed (5,028, 21.04%). In the lncRNAs derived from AS events, 269, 159, 136, 296, and 2 lncRNAs were generated from AA, AD, ES, IR, and MX events, respectively (Fig. 2d).
Predicting the target genes and calculating differentially expression of lncRNAs To investigate the potential functional roles of identi ed rice lncRNAs, we predicted the cis target genes by search for mRNA within 100 Kb upstream and downstream sites of these lncRNAs. Finally, total 2,735 mRNAs associated with 2,383 genes may be regulated by identi ed lncRNAs (Additional le 11: Table   S7). Compared with CK group, the differentially expression lncRNAs (DELs) of rice infected with RBSDV (RBSDV vs CK) and RSV (RSV vs CK) were identi ed based on fold change ≥ 2 and P-value < 0.05 by performing R package DEseq2. There were 344 DELs including 172 up-regulated and 172 down-regulated lncRNAs identi ed in RBSDV vs CK (Fig. 3a), 176 DELs including 61 up-regulated and 115 down-regulated lncRNAs identi ed in RSV vs CK ( Fig. 3d; Additional le 12: Table S8). To further investigate the possible function of DELs in different groups, we performed KEGG pathway and GO term enrichment analysis of targets of DELs by using KOBAS and topGO, respectively (Additional le 13: Table S9). In RBSDV vs CK group, the DELs were associated with 78 KEGG pathways, including 'Carbon metabolism', 'Plant hormone signal transduction', and 'Phenylpropanoid biosynthesis' (Fig. 3b) and 272 GO terms of biological process (BP), including cellular process, biological regulation, and response to stimulus (Fig. 3c). In RSV vs CK group, we found that there were 8, 7, and 5 targets encoding 'Spliceosome', 'Starch and sucrose metabolism', and 'Plant-pathogen interaction', respectively, in enriched KEGG pathways (Fig. 3e) and for GO terms of biological process, major categories were cellular nitrogen compound metabolic process, cellular aromatic compound metabolic process, and organic cyclic compound metabolic process (Fig.  3f).

Characterization of share DELs
To further investigate the share lncRNAs associated with rice response to RBSDV and RSV, we performed distribution analysis of all DELs in both groups. The venn diagram showed that a total of 294 and 126 speci c DELs of RBSDV vs CK and RSV vs CK, respectively, and 50 share DELs of both groups (Additional le 2: Fig. S2). We used R package TCseq to cluster the share DELs based on expression patterns and 4 clusters were generated (Fig. 4a). There were 7, and 15 up-regulated lncRNAs from cluster2, and cluster4, respectively, and 19 down-regulated lncRNAs from cluster3 in both comparisons. The targets of 22 upregulated lncRNAs from cluster2, and 4, were subjected to KEGG pathway and GO term enrichment analysis and multiple KEGG pathways were signi cantly enriched (P < 0.05), such as 'Sesquiterpenoid and triterpenoid biosynthesis', 'Carotenoid biosynthesis', and 'Proteasome'. Meanwhile, several GO terms associated with biosynthetic process were signi cantly enriched (P < 0.05), such as carbohydrate biosynthetic process, polysaccharide biosynthetic process, and glucan biosynthetic process (Fig. 4b).
Verifying accuracy of the RNA-seq by qRT-PCR To verify the RNA-seq data, total 10 co-upregulated lncRNAs from cluster2 and 4 were randomly selected to perform qRT-PCR. The result showed that there were 8 signi cantly up-regulated lncRNAs in rice infected with RBSDV and RSV compared with CK (Fig. 5). The expression tendencies of most of selected lncRNAs were consistent with the sequencing results, indicating that the RNA-seq data was accuracy and the function of these lncRNAs in rice response to RBSDV and RSV were supported.

Construction of regulatory network
To globally identify transcription factors (TFs) in rice, we uploaded sequences of all identi ed transcripts into Plant Transcription Factor Database and nally obtained 1,535 transcripts encoding 55 TFs (Additional le 14: Table S10). In these TFs, bZIP and WRKY were related to most transcripts (112), followed by bHLH (107) (Additional le 3: Fig. S3). To understand the regulatory patterns between TF, lncRNA, and targets, the directed networks were constructed by using R package GENIE3. We rst predicted 17,885 transcripts including 15,969 mRNA, and 1,916 lncRNAs targeted by TFs. Subsequently, we calculated the regulated relationship between the 15,969 mRNAs (target genes) and 1,916 lncRNAs (regulatory factors) by using GENIE3. Finally, total 1,145 co-upregulated transcripts including 60 genes encoding TFs (23 TFs) (regulatory factors), 21 lncRNAs (regulatory factors or targets), and 1,064 mRNAs (targets) were generated (Additional le 15: Table S11). In these transcripts, each lncRNA regulated by corresponding TF can regulate mRNA targeted by the TF, suggesting that the TF may regulated the mRNA via the corresponding lncRNA. The 1,064 mRNAs were subjected to KEGG pathway analysis and 10 pathways were signi cantly enriched (P < 0.05) (Additional le 4: Fig. S4; Additional le 16: Table S12). There were 17, 16, and 12 mRNAs associated with 'Plant hormone signal transduction', 'Phenylpropanoid biosynthesis', and 'Plant-pathogen interaction' (Additional le 4: Fig. S4; Additional le 16: Table S12). We used Cytoscape software to construct directed regulatory network between 60 TFs, 21 lncRNAs, and signi cantly enriched pathways associated (Fig. 6a). In the network, lncRNAs (OS_LNC0494, OS_LNC1144, OS_LNC1812) regulated by Os01t0730700-01 encoding WRKY14 may target ABA responsive element binding factor ABF (Os01t0867300-01) involved in ABA signaling pathway.
Performing qRT-PCR showed that the expression levels of the three lncRNAs were signi cantly upregulated in RBSDV and RSV infected samples (Fig. 6b). We measured the expression of WRKY14, ABF, and the three lncRNAs in rice treated with ABA for 12h and 24h and found that WRKY14, and ABF were signi cantly increased compared with control ( Fig. 6c). In addition, OS_LNC1812 were also up-regulated in treated samples compared with control (Fig. 6c). The result implied that OS_LNC1812 regulated by WRKY14 may be associated with rice response to RBSDV and RSV infection via ABA signaling pathway.

Discussion
The lncRNAs were identi ed in multiple plants, such as Arabidopsis [36], maize [37], rice [38], cotton [39], wheat and Medicago [40]. However, these reports mainly demonstrated the functional characteristics in plant reproductive development and vegetative growth. There were limited data available for the roles of lncRNA in plant response disease and biotic stress. In rice, lncRNAs associated with defense against RBSDV infection were globally identi ed [15,41]. To data, there were no information on lncRNA that respond to RSV infection in rice. This work aimed to identify lncRNAs that respond to RBSDV and RSV by performed strand-speci c paired-end deep sequencing of rice leaf samples infected with the two pathogens and a total of 1,925 lncRNAs were identi ed. The length and expression of lncRNAs were shorter and lower than protein-coding mRNAs in plants [42]. Most lncRNAs were two-exon transcripts, and analysis of TPM distribution also implied the relatively lower expression of lncRNAs compared with mRNAs. The lncRNAs were more evenly distribution across all chromosomes than mRNAs in plants [24]. In this work, we observed most mRNAs were located in the chromosome "arms" of most rice chromosomes and lncRNAs were evenly distributed on 12 rice chromosomes. These characteristics of lncRNAs obtained in this work were consistent with previously reports in rice or other organisms.
Alternative splicing (AS), an important regulatory mechanism of biological process, can generate more than one transcripts including mRNA, and lncRNA [24,43] in eukaryotes by using different splice sites. There were ve fundamental groups categorized from AS events, including intron retention (IR), exon skipping (ES), alternative donor (AD), alternative acceptor (AA) and mutually exclusive exons (MX) [44]. The previous reports had estimated that higher 60% of multiple exon-containing genes were generated from AS events in Arabidopsis [45]. In this work, we detected 23,898 AS events associated 724 lncRNAs and 4,424 mRNAs. Compared other four groups, IR events represented the largest proportion of all AS events (27.14%) and the number of lncRNAs with IR events was estimated to be 44.88% of all ASgenerated lncRNAs. The similar proportion of IR events in this work was observed in other studies in plants [46,47].
LncRNAs mediated gene expression and protein synthesis by multiple ways, such as transcription, chromatin modi cation, and post-transcriptional processing [48,49]. Therefore, the critical ways to elucidate the function of lncRNAs is predication of LncRNA target genes and total 5,128 cis-target mRNA associated with 4,379 genes were identi ed. Difference analysis showed that 344 and 176 DELs in RBSDV vs CK and RSV vs CK, respectively, were identi ed. The targets of DELs from the two comparisons were subjected to KEGG pathway and GO term enrichment analysis, respectively and multiple resistancerelated pathways were signi cantly enriched, such as 'Plant-pathogen interaction', and 'Plant hormone signal transduction' [50,51]. Distribution analysis of all DELs indicated that total 50 share DELs between both groups. All share DELs were divided into four clusters by using R package TCseq and 22 DELs of cluster2, and 4 were co-upregulated. The targets of the 22 DELs were associated with many pathways, such as 'Sesquiterpenoid and triterpenoid biosynthesis', 'Carotenoid biosynthesis', and 'Proteasome'.
Asparagine-rich protein (NRP) of Nicotiana benthamiana (NBnrp) can regulate sesquiterpenoid biosynthesis to participate in defense responses induced by PevD1, a fungal effector of Verticillium dahlia [52]. The proteasome targeted by multiple bacterial effectors is an important component of systemic acquired resistance (SAR) and pathogen-associated molecular pattern-triggered immunity [53]. In addition, we used qRT-PCR to measure the expression of the 22 share DELs and found that the expression level of all DELs were increased in RBSDV and RSV infected samples compared with control. The above results suggested that these share up-regulated DELs may be involved in rice response to RBSDV and RSV infection. To further investigate the mechanism, we constructed the gene regulatory network (GRN) to explore the regulatory patterns of 1,535 genes encoding TF (regulatory factors), 22 common up-regulated DELs (regulatory factors or targets) and mRNA (target genes) during RBSDV and RSV infection in rice by using R package GENIE3. Finally, total 21 share up-regulated lncRNAs regulated by 60 genes encoding TF can regulate 1,064 co-up regulated mRNAs targeted by the 60 TF genes. Functional enrichment analysis of the 1,064 mRNAs showed that multiple defense-related pathways were signi cantly enriched, such as'Plant hormone signal transduction', 'Phenylpropanoid biosynthesis', and 'Plant-pathogen interaction'.
Phytohormones, such as ABA, have been shown to response against microbial pathogens in plant [54]. In Arabidopsis, the transcription factor WRKY33 can increase the transcription of NCED3 and NCED5 to promote ABA biosynthesis after necrotrophic pathogen Botrytis cinereal infection [55]. ABA increases susceptibility to RBSDV infection by inhibiting the jasmonate pathway and regulating reactive oxygen species levels in rice [13]. Construction of GRN showed that OS_LNC1812 regulated by Os01t0730700-01 encoding WRKY14 can target ABF (Os01t0867300-01) maintaining new homeostatic levels by creating a feedback control loop [56]. Performing qRT-PCR showed that the expression levels of WRKY, ABF and OS_LNC1812 were increased in rice after ABA treatment and infected with RBADV and RSV, respectively.
These results implied that ABA, RBSDV, and RSV may elevate expression ABF by WRKY14-induced OS_LNC1812, and the high level of ABF can decrease ABA synthesis to increase rice resistance. Although we described OS_LNC1812 may be associated with ABA signaling pathway, the molecular mechanism behind this of response to pathogens infection was still unclear. These concerns give rise to possible clues exploring the defense mechanism of the OS_LNC1812-modulated ABA pathway. The epigenetic regulation of this lncRNA was needed to further investigate at the molecular level.

Conclusion
Our study systematically identi ed lncRNA expression pro les of rice infected with RBSDV and RSV, respectively. We obtained a total of 1,925 lncRNAs including 1,112 genic lncRNAs, and 813 intergenic lncRNAs from RNA-seq data. The lncRNAs were shorter, lower expression, and more evenly distribution than mRNAs. Analyzing alternative splicing events indicated that total 724 lncRNAs and 4,424 mRNAs were generated from 23,898 AS events. Differentially expression analysis showed that there were 344, and 176 DELs generated from RBSDV vs CK and RSV vs CK, respectively. The targets of DELs were involved with multiple defense-related pathways, such as 'Plant hormone signal transduction', 'Plantpathogen interaction', and 'Phenylpropanoid biosynthesis'. In addition, distribution analysis of DELs and construction of gene regulatory network (GRN) indicated that total 21co-upregulated lncRNAs regulated by 60 TFs may target 1,064 mRNAs regulated by the 60TFs. In the GRN, lncRNAs (OS_LNC0494, OS_LNC1144, OS_LNC1812) regulated by Os01t0730700-01 encoding WRKY14 (Os01t0730700-01) may target ABF (Os01t0867300-01) and the three lncRNAs, Os01t0730700-01, and Os01t0867300-01 were selected for further research. Our study provided genome-wide lncRNA expression pro les and novel potential modulation of lncRNAs involved in plant disease resistance via ABA pathway in rice.

Plant material, virus inoculation and RNA isolation
The rice cultivar LinDao10 was used in this study and the seedlings were grown inside a greenhouse at 28 •C and 70%-80% humidity with a photoperiod of 16h light and 8 dark. The rice seedlings at 3-leaf-tage were inoculated with RBSDV and RSV by using viruliferous small brown planthoppers (SBPH, Laodelphax striatellus). The viruliferous SBPHs were feed using assayed seedlings for two days and were then removed. We used rice seedlings inoculated with non-viruliferous SBPHs as controls (CK) in this work.
The leaves with disease phenotype in experiment groups and health leaves in control were harvested at 30 days after the inoculation. Subsequently, all collected leaves were quickly frozen in liquid nitrogen, and then stored at -80 •C separately until use. We used the mixed leaf tissues collected from 10 plants of each group as a biological replicate, and each treatment had three biological replicates. AG RNAex Pro Reagent AG21101 (Accurate Biotechnology (Hunan) Co., Ltd) was used to extract the total RNA of individual samples, and the speci c primers for RBSDV and RSV were designed to con rm successfully virus infection in each virus-inoculated sample by using RT-PCR (Table S1).

Preparation RNA library and sequencing
The transcriptome library construction and Illumina deep sequencing were performed by OE Biotech Co., Ltd. (Shanghai, China). Epicentre Ribo-Zero TM Kit (Epicentre, Madison, WI, USA) was used to remove Ribosomal RNA. Subsequently, we used fragmented RNAs as templates to amplify cDNA with random hexamers. DNA polymerase I and RNase H were used to synthesize second-strand cDNA. Illumina Hiseq™ 4000 sequencing system (Illumina, China) was performed to sequence the library and reads with 150 bp paired-ends were generated. We deposited the raw data in Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) and the accession number was PRJNA699905.
LncRNA identi cation and target gene prediction Trimmomatic software (v. 0.39) was performed to lter the low-quality bases containing adapter or poly-N and high-quality clean reads were generated [57]. The spliced read aligner HISAT (v2.1.0) [58] was used to align all clean reads to the rice reference genome (IRGSP1.0) [59] and StringTie software (v2.1.3) was used to assemble all transcripts [60]. Gffcompare (v. 0.11.2) was used to compare assembled transcript isoforms to potato genome annotation [61]. Salmon (v 0.14.1) software was used to calculate the expression level of all assembled transcripts based on calculating Transcripts Per Million (TPM) values [62]. Tximport (v 1.18.0) was used to import transcript-level abundance and summarize into expression matrices for downstream analysis. All assembled transcripts were used to identify lncRNAs by using FEELnc software (v. 0.1.1) [63], CPC [64], and Swiss-Prot database [65]. First, FEELnc software was used to lter short transcripts (< 200 bp), single exon transcripts, transcripts with ORF greater than 100 amino acids and compute a coding potential score. Second, the Swiss-Prot database and CPC software were used to calculate the coding potential and transcripts with potential coding capabilities were removed. Third, the expression level of transcripts with TPM < 1 in all samples were ltered. Finally, total 1,925 transcripts were identi ed as lncRNAs. We de ned the coding genes within 100 kb upstream and downstream of lncRNAs as the cis-targets. To further understand the biological process of lncRNAs, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis were performed by using KOBAS software 3.0 [66] and topGO software (v 2.42.0) [67], respectively.

Differentially expression analysis
The transcripts with TPM > 1 at least one sample were considered to be expressed. The expressed transcripts including mRNAs and lncRNAs were used to calculate difference by performing R package DEseq 2 (v1.30.0) with P-value ≤ 0.05 and log 2 (fold change) > 1 [69]. We performed R packages circlize(), ggplot(), and TCseq to draw distributions of lncRNAs and mRNAs on rice chromosomes, bubble charts, heatmaps, and expression patterns maps, respectively. TBtools software (v1.051) was used to plot venn gram [70]. Cytoscape software (v3.7.2) was used to construct interactive network between lncRNAs, TFs, and pathways associated with targets [71]. R package GENIE3 (v1.12.0) was used to construct gene regulatory network (GRN) [72]. Brie y, the mRNAs, TFs, and lncRNAs with TPM > 1 at least one sample were selected to construct GRN based on expression patterns. Subsequently, the connectivity between regulatory factor (TF or lncRNA) and targets (lncRNA or mRNA) was calculated by GENIE3 software. The interacted pair with connectivity > 0.3 were extracted [73].

Real-time quantitative PCR analysis
Real-time RT-PCR was performed using LightCycler96 (Switzerland). We used Evo M-MLV RT Kit with gDNA Clean for qPCR AG11705 (Accurate Biotechnology (Hunan) Co., Ltd) to reverse-transcribe the total RNA into a single-stranded RNA. We designed the speci c primers by using Primer-BLAST in NCBI (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). We selected ACTIN2 to be an internal reference gene.
We performed real-time PCR analysis following the manufacturer's instructions of SYBR PrimeScript TM RT-PCR Kit (Takara, China) with three biological replicates and each biological replicate contained three technical replicates. Relative expression levels of detected mRNA or lncRNA were calculated using the 2 -ΔΔCt method [74]. All primers used in this work were listed in table S1.

ABA treatment
We used ABA treatments in the same method as described by Liu et al [75]. Briefly, we used 100 µM ABA solutions to spray rice seedlings at the three-to four-leaf stage and rice sprayed with dimethylsufoxide (DMSO) were used as the controls. After treatment for 12h and 24h, the rice samples were collected. The 10 plants were used of each treatment as a biological replicate, and there were three biological replicates for each treatment.

Statistical analysis
The data were represented by mean±standard deviation (SD). For analyzing the differences between groups, the one-way analysis of variance (ANOVA) and Tukey's test were used. P-value ≤ 0.05 was considered as statistically signi cant.      Validation of 10 co-upregulated expressed lncRNAs from cluster by qRT-PCR. The values were represented by means of the gene expression levels ± standard deviations (SD) (n = 3 biological replicates). *P < 0.05, ***P < 0.001 represented statistical signi cance.

Figure 6
Construction of gene regulatory network using transcription factors, mRNAs, and lncRNAs. (A) The 21 coupregulated lncRNAs regulated by 60 genes encoding TFs targeted mRNAs-related pathways. (B) QRT-PCR was used to measure the expression levels of lncRNAs. (C) Measuring the expression levels of WRKY14, ABF, and the three lncRNAs in rice treated by ABA. The values were represented by means of the gene expression levels ± standard deviations (SD) (n = 3 biological replicates). * P < 0.05, **P <0.01, ***P < 0.001 represented statistical signi cance.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.