Systematic Characterization of Long Non-Coding RNAs and Their Responses to Drought Stress in Dongxiang Wild Rice

Qi Weidong, Chen Hongping, Yang Zuozhen, Hu Biaolin, Luo Xiangdong, Ai Bing, Luo Yuan, Huang Yu, Xie Jiankun, Zhang Fantao

Research Paper

Systematic Characterization of Long Non-Coding RNAs and Their Responses to Drought Stress in Dongxiang Wild Rice

Qi Weidong1, #, Chen Hongping2, #, Yang Zuozhen3, #, Hu Biaolin2, Luo Xiangdong1, Ai Bing1, Luo Yuan1, Huang Yu1, Xie Jiankun1, Zhang Fantao1

(College of Life SciencesChina; Rice Research Institute, Jiangxi Academy of Agricultural Sciences, Nanchang 330200, China; Cipher Gene Limited Liability Company, Beijing 100080, China; These authors contributed equally to this work)

Long non-coding RNAs (lncRNAs) play important roles in response to various biotic and abiotic stresses. So far, systematic identification and characterization of lncRNAs have been reported in a few model plant species and major crops, but their roles in abiotic stress response have not yet been reported in common wild rice (). Dongxiang wild rice (DXWR) possesses a high degree of drought resistance and has been well recognized as a precious genetic resource for drought resistant rice breeding. We presented the reference catalog of 1 655 novel lncRNA transcripts in DXWR using strand-specific RNA sequencing and bioinformatics approaches. Meanwhile, a total of 1 092 lncRNAs were determined as differentially expressed lncRNAs under drought stress. Quantitative real-time PCR results exhibited a high concordance with RNA sequencing data, which con?rmed that the expression patterns of lncRNAs based on RNA sequencing were highly reliable. Furthermore, 8 711 transcripts were predicted as target genes of the differentially expressed lncRNAs. Functional annotation analysis based on Gene Ontology and Kyoto Encyclopedia of Genes and Genomes databases showed that the target genes were signi?cantly enriched in cellular and metabolic processes, cell part, binding and plant hormone signal transduction, as well as many other terms related to abiotic stress resistance. These results expanded our understanding of lncRNA biology and provided candidate regulators for genetic improvement of drought resistance in rice cultivars.

Dongxiang wild rice;drought stress; long non-coding RNA; systematic characterization

Rice () is a major cereal crop in many parts of the world, serving as the staple food for more than half of the world’s population (Yang et al, 2015). Due to global warming and increasing scarcity of water resource, rice production has been affected seriously by drought stress(Shanker et al, 2014; Duan et al, 2019). Therefore, breeding rice cultivars with drought resistance is a promising strategy to meet the world’s increasing demand for food (Luo, 2010). However, it is particularly challenging since drought resistanceis a highly complex trait mediated by an elaborate network of genetic pathways(Anupama et al, 2018).

Over the last 20 years, although many related genes have been cloned and functional characterized in different plant species (Joshi et al, 2016; Moran et al, 2017), the mechanisms underlying drought resistance remain elusive. The discovery of non-coding RNAs (ncRNAs) provides a new insight into genome regulation (Heo et al, 2013). Based on the length, ncRNAs can be sub-divided into small ncRNAs (less than 200 nucleotides) and long ncRNAs (lncRNAs, more than 200 nucleotides) (Cao, 2014). Till now, the roles of small ncRNAs (particularly miRNAs) have been extensively studied in plants because of their function as gene regulators during plant growth and development, as well as in response to various biotic and abiotic stresses (Sun et al, 2014; Li et al, 2017; Liuet al, 2018). Following the advance of sequencing technology, although systematic identification of lncRNAs has been carried out in a few major crops and model plant species, there are still a myriad of lncRNAs remain unidentified in other plant species or closely related plant species (Zaynab et al, 2018).

As the progenitor of cultivated rice, common wild rice () is a valuable plant germplasm resource that has many specific characters (Menguer et al, 2017). Dongxiang wild rice (DXWR) is a common wild rice that originated from Dongxiang county (28o14N, 116o36E), Jiangxi Province, China, which is the northernmost region in the world where.has been found to date (Zhang et al, 2016; Liang et al, 2018). Previous studies have revealed that DXWR can survive under extreme drought stress condition, suggesting that it is a precious plant germplasm resource for genetic improvement of drought resistance in cultivated rice (Zhang et al, 2006). However, the molecular mechanisms underlying drought resistance of DXWR remain unclear.

In this study, we used strand-specific RNA sequencing and bioinformatics approaches to identify and characterize lncRNAs in DXWR under drought stress and normal conditions. The main objectives were to determine how many and which lncRNAs were differentially expressed in DXWR under drought stress, which genes were predicted to be targeted by the differentially expressed lncRNAs, and which biological processes and pathways were significantly enriched for the predicted target genes. This may lay the foundation for future functional and evolutionary studies on lncRNAs for common wild rice in response to drought stress.

Materials and Methods

Plant materials and drought stress treatment

The experiments were conducted in Jiangxi Normal University (28o40N, 115o55E, Nanchang, Jiangxi Province, China) in the summer of 2017. Drought stress treatment was performed following the methods of Zhang et al (2016). The seeds of DXWR from Shuitaoshuxia population were immersed in distilled water in the dark for germination, and the uniformly germinated seeds were then grown in a plastic pot with IRRI (International Rice Research Institute) nutrient solution (1.25 mmol/L NH4NO3, 0.3 mmol/L KH2PO4, 0.35 mmol/L K2SO4, 1 mmol/L CaCl2·2H2O, 1 mmol/L MgSO4·7H2O, 0.5 mmol/L Na2SiO3, 20 μmol/L NaFeEDTA, 20 μmol/L H3BO3, 9 μmol/L MnCl2·4H2O, 0.32 μmol/L CuSO4·5H2O, 0.77 μmol/L ZnSO4·7H2O, and 0.39 μmol/L Na2MoO4·2H2O, pH 5.5). The plastic pot with seeds was then transferred into a plant growth chamber, and the parameters were set as follows: the day/night temperature was 30 oC/ 26 oC (14 h day/10 h night), and the relative humidity was 70%. The nutrient solution in the plastic pot was renewed every 3 d. For drought stress treatment, when the plants were at four-leaf stage, half of them were exposed in the air of the growth chamber under light for 12 h, and the leaves of plants were obviously wrinkled. The untreated plants at the same time point were used as control. For each sample, six plants were harvested and pooled in order to minimize the effect of transcriptome unevenness among plants. The collection was repeated three times as biological replicates and immediately frozen in liquid nitrogen and then held at -80 oC.

RNA isolation, library construction and sequencing

Total RNA isolation, library construction and strand- specific RNA sequencing were performed by the Annoroad Gene Technology Corporation (Beijing, China). Total RNA was isolated using TRIzol reagent following the manufacturer’s protocol (Qiagen, Valencia, USA). RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using a Kaiao Photometer Spectrophotometer K5500 (Kaiao, Beijing, China). RNA integrity and concentration were measured with the RNA Nano 6000 Assay Kit in an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, USA). Approximately 3 μg of total RNA for each sample was used to construct the library. The ribosomal RNA (rRNA) was removed by Ribo-ZeroTMGold Kit (Illumina, San Diego, USA). NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ispawich, USA) was used for library preparation and all of the procedures and standards were performed according to the manufacturer’s protocol. After quality control, the library preparations were sequenced on an Illumina HiSeq XTen platform (Illumina Inc., San Diego, USA) and 300 bp paired- end reads were generated.

Raw data filtering

Before further analysis, raw data filtering was needed to decrease data noise. The steps were as follows: reads containing the adapter sequences were removed; low-quality reads (reads containing more than 15% bases with≤ 19) were removed; reads containing more than 5% unknown nucleotides were removed. Then, in consideration of rRNA pollution interference in the analysis, the remaining reads were mapped to the rRNA reference sequence using the short-read alignment software SOAP2 (http://soap.genomics.org.cn/) to remove the rRNA reads. Finally, the remaining reads were named clean reads and used for further analyses.

Quality control, data processing and differential expression analysis

Quality control and read statistics were determined by FastQC (Kroll et al, 2014). The clean reads were mapped to the.ssp.cv. Nipponbare reference genome using HiSAT2 (Kim et al, 2015). All lncRNAs and mRNAs detected in the RNA sequencing were obtained by StringTie (Pertea et al, 2015). To reduce the false positive rates, the candidate lncRNAs had to meet the following requirements: transcripts should contain more than one exon and be no shorter than 200 bp; transcripts should have more than five reads of coverage; transcripts belong to coding-genes, pseudogenes, rRNA, tRNA, snoRNA and snRNA were removed. Protein coding potency of transcripts was analyzed by four programs: Coding- Non-Coding Index (CNCI) (Sun et al, 2013), Coding Potential Calculator (CPC) (Kong et al, 2007), Protein Families Database (PFAM) (Punta et al, 2012), and Coding-Potential Assessment Tool (CPAT) (Wang et al, 2013). LncRNA and mRNA sequence reads were normalized to FPKM (fragments per kilobase per million mapped fragments) (Trapnell et al, 2010). Differential expression analysis was performed using the DEGseq R package (Wang et al, 2010). An absolute value of log2(fold change) ≥ 1 and< 0.05 were set as the threshold for significant differential expression.

Quantitative real-time PCR (qRT-PCR)

The reverse transcription of total RNA was performed according to the instruction of the M-MLV Reverse Transcription Kit (Takara, Dalian, China). The qRT-PCR was conducted using SYBR Premix ExII (Takara, Dalian, China) and ABI 7500 Real-time PCR System (Applied Biosystems, Carlsbad, USA). The primers were synthesized by Sangon Biotechnology (Shanghai, China) and listed in Supplemental Table 1. The housekeeping genewas used as the internal control. The relative expression levels were calculated using 2-ΔΔC(t)method as described by Jain et al (2006).

Target gene prediction and enrichment analysis

The lncRNAs act on neighboring target genes, which is known as therole of lncRNAs (Zaynab et al, 2018). The transcripts at 50 kb upstream and downstream were searched to predict putative target genes of the differentially expressed lncRNAs, followed by analyzing the functions of these transcripts to annotate the lncRNAs. Gene Ontology (GO) enrichment analysis of lncRNA target genes was performed using the GOseq R package (Young et al, 2010). In addition, the KOBAS software (Mao et al, 2005) was used to test the statistical enrichment of target genes in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GO terms and KEGG pathways with< 0.05 were considered significantly enriched.

Results

Basic data statistics

In total, 116 497 948 and 103 360 478 raw reads were obtained in the control (DY-CK) and drought stress (DY-D) libraries, respectively (Table 1). After trimming adapters and filtering out low quality reads, 92 556 246 (DY-CK) and 83 242 176 (DY-D) clean reads were remained and used for further analysis (Table 1). To estimate the quality of sequencing data, we used FastQC programe to calculate the quality score () of each base pair in the reads. The meanvalues of both DY-CK and DY-D libraries were 40, showing that the sequencing data were highly reliable (Supplemental Fig. 1). A total of 87 958 100 and 78 406 778 clean reads in the DY-CK and DY-D libraries, respectively, were mapped to the reference genome (Table 1). According to the percent of reads mapped to genomic regions, we found approximately 66.01% reads were mapped to exonic region in the DY-CK librarywhile an increase was observed in the DY-D library (71.79%) (Table 1). In contrast, the proportions of intronic and intergenic reads were lower in the DY-Dlibrary compared with those in the DY-CK library (Table 1).

Table 1. RNA sequencing data for two samples.

Values in the parenthesis are the ratio to raw read number;Values in the parenthesis are the ratio to clean read number.

Identi?cation and characterization of lncRNAs in DXWR

In total, 1 655 novel lncRNA transcripts with high reliability corresponding to 1 461 lncRNA genes were identified and confirmed by the CNCI, CPC, PFAM and CPAT programes (Fig. 1 and Supplemental Table 2). Meanwhile, the majority (74.80%) of lncRNAs were located in the intergenic regions (SupplementalTable 2). Similar to mRNAs, lncRNAs were widespread in all chromosomes. The number of lncRNAs varied among different chromosomes with the maximum number (291) in chromosome 1 and the minimum number (94) in chromosome 7 (Supplemental Table 2). To study the basic features of lncRNAs, we compared the lncRNAs with mRNAs and found that lncRNAs were shorter in transcript length. The length of lncRNAs varied from 200 to 4 037 nucleotides with an average length of 527 nucleotides (Supplemental Fig. 2 and SupplementalTable 3). The majority of lncRNAs (90.51%) were short-length lncRNAs (< 1 000 nucleotides), and 7.79% were medium-length lncRNAs (1 000–2 000 nucleotides), and only 1.69% were long-length lncRNAs (> 2 000 nucleotides) (Supplemental Fig. 2 and SupplementalTable 3). In contrast, medium-/long-length mRNAs (≥ 1 000 nucleotides) accounted for 67.03% of the total number of mRNAs (Supplemental Fig. 2 and SupplementalTable 4). Comparison of the exon features of lncRNAs and mRNAs showed that lncRNAs contained 2.3 exons on average per transcript, which was fewer than that of mRNAs (4.5 exons on average) (Supplemental Tables 3 and 4). Most (78.85%) lncRNAs contained two exons (Supplemental Fig. 2 and SupplementalTable 3), while two-exon mRNAs accounted for only 14.65% of the total mRNAs (Supplemental Fig. 2 and SupplementalTable 4), and the number was apparently lower than that of lncRNAs. Moreover, we analyzed the average expression levels of lncRNAs and mRNAs, and boxplots illustrated that the expression levels of lncRNAs were significantly lower than those of mRNAs in both DY-CK and DY-D libraries (Supplemental Fig. 3).

Fig. 1. Venn diagram of non-coding transcripts from the Coding-Non-Coding Index (CNCI), the Coding Potential Calculator (CPC), the Protein Families Database (PFAM), and the Coding-Potential Assessment Tool (CPAT).

Differential expression of lncRNAs in response to drought stress

Based on the sequencing results, the lncRNAs with log2(fold change) ≥ 1 and< 0.05 were designated as signi?cantly up-regulated, whereas the lncRNAs with log2(fold change) ≤ -1 and< 0.05 were designated as signi?cantly down-regulated. The statistical results indicated that a total of 1 092 lncRNAs were differently expressed in response to drought stress, where 542 lncRNAs signi?cantly up-regulated and 550 lncRNAs were down-regulated (Supplemental Table 5). Meanwhile, 420 lncRNAs were co-expressed in the two libraries, while 328 and 344 lncRNAs were specific expressed in the DY-CK and DY-D libraries, respectively (Supplemental Table 5). The most signi?cantly up-regulated lncRNA was MSTRG69391, followed by MSTRG41712 and MSTRG68635. Conversely, the most signi?cantly down-regulated lncRNA wasMSTRG65848, followed by MSTRG27834 and MSTRG46301 (Supplemental Table 5).

To confirm the expression patterns of lncRNAs obtained by RNA sequencing, we used qRT-PCR method to analyze the expression of 15 randomly selected lncRNAs with differently expressed levels. According to the sequencing results, these lncRNAs included seven down-regulated (MSTRG25356, MSTRG29971, MSTRG42925, MSTRG45354, MSTRG54123, MSTRG55730 and MSTRG69646) and eight up-regulated (MSTRG5459, MSTRG14660, MSTRG21367, MSTRG22525, MSTRG31874, MSTRG43289, MSTRG46940 and MSTRG62647) lncRNAs. As a result, the expression patterns of lncRNAs detected by qRT-PCR were highly consistent with those based on RNA sequencing (Fig. 2), which indicated that the expression patterns of lncRNAs based on RNA sequencing were highly reliable.

Prediction of target genes for drought-responsive lncRNAs

A total of 8 711 transcripts were predicted as the target genes of the drought-responsive lncRNAs(SupplementalTable 6). As for the number of target genes per lncRNA, MSTRG42334 had the maximum number (34) of target genes, whereas no target genes were predicted for MSTRG17879, MSTRG17882, MSTRG21009, MSTRG21492, MSTRG40669, MSTRG51860, MSTRG51861, MSTRG51865, MSTRG55763, MSTRG56497, MSTRG70414, MSTRG72630, MSTRG72631, MSTRG72672 and MSTRG72686 (Supplemental Table 7). Meanwhile, 5 010 transcripts were significantly differently expressed under drought stress, where 2 282 and 2 728 transcripts were significantly up- and down-regulated, respectively (Supplemental Table 6).Many of these target genes were predicted to encode a variety of proteins, such as transcription factor (TF) proteins, transposon and retrotransposon proteins, growth-regulating factor proteins, RNA-binding proteins, mitochondrial phosphate carrier proteins and calmodulin-binding heat-shock proteins. These results implied that the differentially expressed lncRNAs can play vital regulatory functions in a wide range of biological processes to cope with drought stress condition.

GO enrichment and KEGG pathway analyses for potential target genes

To delineate the biological roles, the predicted target genes of the drought-responsive lncRNAs were sorted into GO term categories. The results illustrated that the target genes can be clustered into 1 088, 2 786 and 6 840 GO terms in cellular component (Supplemental Table 8), molecular function (Supplemental Table 9) and biological process (Supplemental Table 10) categories, respectively. Within the biological process category, cellular and metabolic processes were the most highly represented groups, which suggested that extensive metabolic activities were taking place in the drought treated plants (Fig. 3). Within the cellular component category, transcripts that corresponded to cell part, organelle and membrane were the most abundant, and binding, catalytic and transporter were the most abundant groups within the molecular function category (Fig. 3). Meanwhile, the target genes were also found to be enriched in other stress- related terms, such as response to stimulus, biological regulation, nucleic acid binding transcription factor, and molecular function regulator (Fig. 3). Furthermore, functional hierarchy of the target transcripts in the KEGG orthology system can be classi?ed into five categories, namely, environmental information processing, genetic information processing, metabolism, cellular processes and organismal systems (Fig. 4). The statistical results demonstrated that most of the target genes were related to metabolism, followed by genetic information processing and organismal systems (Fig. 4). Moreover, a KEGG reference pathway analysis was carried out. As shown in Supplemental Table 11, the significantly enriched pathways (< 0.05) were plant hormone signal transduction(= 0.0008), ABC transporters (= 0.0152), flavonoid biosynthesis (= 0.0192), one carbon pool by folate (= 0.0218), benzoxazinoid biosynthesis (= 0.0366), and tropane, piperidine and pyridine alkaloid biosynthesis (= 0.0442).

Fig. 2. Comparison of expression levels of 15 randomly selected drought responsive lncRNAs using RNA sequencing and qRT-PCR methods.

Fig. 3.Gene ontology (GO) classification for the predicted target genes of drought-responsive lncRNAs.

Fig. 4.Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway assignments for predicted target genes of drought-responsive lncRNAs.

Discussion

LncRNAs, as a type of non-coding RNA molecules, are increasingly being unraveled in plants, animals, and fungi (Zhang et al, 2017). However, few studies have reported the detailed information of lncRNAs in common wild rice. In the present study, we identi?ed and characterized lncRNAsin DXWR on a genome- wide scale. Consistent with previous studies (Li et al, 2014;Zhu et al, 2015), most of the lncRNAs identi?ed in the present study were located in the intergenic regions, indicating that the type of lncRNA was long intergenic non-coding RNA. Comparative analysis of lncRNAs and mRNAs revealed the characteristics of them were similar to those of previous studies. For example, previous studies showed that lncRNAs are shorter and have lower expression levels than mRNAs (Pauli et al, 2012; Li et al, 2014). In the present study, the average length of mRNAs was 1.81-fold longer than that of lncRNAs, and the expression levels of lncRNAs were significantly lower than those of mRNAs.

Till now, lncRNAs participating in drought-responsive regulation have been studied in maize (Zhang et al, 2014), cotton (Lu et al, 2016),(Qin et al, 2017), cassava (Li et al, 2017), wheat (Cagirici et al, 2017), cultivated rice (Yuan et al, 2018), etc., implying that lncRNAs as ubiquitous regulators are involved in responding to drought stress in various kinds of plant species. In this study, a total of 1 092 differentially expressed lncRNAs were identified as drought stress- responsive lncRNAs in DXWR.Meanwhile, 8 711 transcripts were predicted as candidate target genes for the drought-responsive lncRNAs, including 2 282 and 2 728 significantly up- and down-regulated transcripts under drought stress. It is noteworthy that some differentially expressed target genes were previously reported to be involved in drought resistance, such as() (Xiong and Yang, 2003),() (Seo et al, 2011),() (Ye et al, 2011),() (Wang et al, 2014),() (Jisha et al, 2015),() (Lee et al, 2015)and() (Ahmad et al, 2016), which have not been connected to lncRNAs yet.

Meanwhile, we also noticed that the other putative target genes were predicted to encode a broad range of proteins and many of them were classi?ed as TFs and functional proteins. Among TF families, the bZIP family is one of the largest and diverse groups (Joshi et al, 2016). In plants, bZIP TFs have been reported to be involved in response to biotic and abiotic stresses (Dr?ge-Laser et al, 2018). Yu et al (2011) revealed that the overexpression of() can enhance heat and drought resistance in rice. In the present study, the expression ofwas significantly up-regulated in DXWR under drought stress. Meanwhile,was predicted to be targeted by lncRNA MSTRG62341, and the expression of MSTRG62341 was significantly down-regulated under drought stress. The regulatory mechanisms of these still need to be expounded in detail. In addition, we found 16 WRKY genes were predicted to be targeted by the drought-responsive lncRNAs. As we all know, the WRKY gene family is another one of the largest TF families in higher plants (Phukan et al, 2016). Till now, many WRKY genes involved in drought stress responses have been identified in different plant species. Ma et al (2017) identifiedin wheat and found that the over-expression ofincan enhancedrought resistance by the induction of stomatal closure. Kiranmai et al (2018) isolatedfrom horsegram and found that the over-expression ofin groundnut leads to a greater resistance to drought stress. An intriguing observation is that, 9 of our predicted 16 WRKY genes were significantly differentially expressed under drought stress, including 8 significantly up-regulatedgenes (,,,,,,and) and 1 significantly down-regulated gene. Meanwhile, according to the previous reports, three of the nine differentially expressed genes are involved in response to abiotic stress, i.e.() (Zhang et al, 2015),() (Song et al, 2009), and() (Song et al, 2010). Certainly, further investigation is required to confirm, if any, the roles of these TF genes and their corresponding lncRNAs involving in the drought resistance of DXWR.

Transposable element (TE) genes can help plants to respond and adapt to various abiotic stresses (Song and Cao, 2017). Yasuda et al (2013) reported that transposoncan induce the expression ofunder salt stress. Meanwhile, the over-expression ofleads to a greater resistance to salt stress than the wild type plants (Xu et al, 2008). Therefore, the transposonhas the potential to improve salt resistance. In this study, 18 TE genes were predicted to be targeted by the drought-responsive lncRNAs. Additionally, the AT-hook motif nuclear localized (AHL) gene family has been reported to play important roles in response to various environmental stimuli (Zhao et al, 2014). Zhou et al (2016) identifiedin rice and found that the over-expression ofenhances drought resistance in rice plantsduring both seedling and panicle developmental stages (Zhou et al, 2016). In this study, five AHL genes (,,,and) were predicted to be targeted by six drought-responsive lncRNAs, including MSTRG30786, MSTRG31067, MSTRG33464, MSTRG45899, MSTRG45897 and MSTRG54431. Among these five AHL genes, the expression levels of two genes (and) were significantly down-regulated and that ofwas significantly up-regulated under drought stress.

Although DXWR has been considered as a precious genetic resource for drought-resistant rice breeding, the molecular mechanisms underlying drought resistance of DXWR are still insufficiently understood (Zhang et al, 2016; Liang et al, 2018). Till now, only 12 drought resistance QTLs have been identified in DXWR (Zhang et al, 2006). Among them,andare the two most significant QTLs for drought-resistance, explaining up to 7% and 14% of the phenotypic variances, respectively (Zhang et al, 2006). In this study, a total of 125 and 65 target genes of drought-responsive lncRNAs were localized within theandintervals, respectively. Based on the annotation of the target genes, we noticed that some of them can be involved in response to biotic and abiotic stresses. The over-expression of zinc-finger protein genefromcan enhance cold and drought resistances in transgenic(Luo et al, 2012) and the over-expression of the zinc-finger protein geneinresults in a greater resistance to drought stress (Yin et al, 2017). These results suggested that the zinc- finger proteins play vital roles in plant responses to a wide spectrum of stress conditions. In the present study, we found that two target genes encoding zinc- finger proteins were located in theinterval, i.e.and. Among them,was significantly up-regulated under drought stress. Meanwhile,was predicted to be targeted by lncRNA MSTRG24141 that was significantly down-regulated under drought stress. Furthermore, miRNAs can involve in the drought resistance of DXWR and a list of target genes for the drought-responsive miRNAs were predicted (Zhang et al, 2016, 2018). Among these target genes, 39 genes were also predicated to be targeted by the drought-responsive lncRNAs, including six significantly up-regulated genes (,,,,and) and eightsignificantly down-regulatedgenes(,,,,,,and) under drought stress. Meanwhile, based on the functional annotation, five differentially expressed target genes can be involved in response to biotic and abiotic stresses, i.e.,,,and. Obviously, these differentially expressed target genes are worthy of further studies.

GO enrichment and KEGG pathway analyses of lncRNAs target genes can help us to understand the functions of lncRNAs more effectively. Our functional prediction according to GO categories showed that the target genes of drought-responsive lncRNAs were significantly enriched in molecular functions including in oxidoreductase and chlorophyllideoxygenase activity. According to Munné-Bosch et al (2009) and Huang et al (2014), lncRNAs can be involved in drought resistance by regulating target genes that control oxidoreductase and chlorophyllideoxygenase activity. In addition, plant hormones and their signal transductionplay important roles in response to various biotic and abiotic stresses (Nawaz et al, 2017). Based on the annotation of the target genes, we found 27, 24, 5, 5 and 5 target genes of drought-responsive lncRNAs were associated with ethylene, auxin, abscisic acid, gibberellin and cytokinin, respectively. Meanwhile, KEGG pathway analysis of the lncRNA target genes indicated that the most significantly enriched pathway was involved in plant hormone signal transduction. Therefore, these results provided an abundant resource of candidate lncRNAs associated with drought resistance and enriched its regulatory network in plants.

AcknowledgementS

This research was supported by the National Natural Science Foundation of China (Grant No. 31660386), the Natural Science Foundation of Jiangxi Province for Distinguished Young Scholars (Grant No. 20171BCB23040), the Foundation of Jiangxi Educational Committee (Grant No. GJJ170193), and the Sponsored Program for Distinguished Young Scholars in Jiangxi Normal University, China.

supplemental data

The following materials are available in the online version of this article at http://www.sciencedirect.com/science/ journal/16726308; http://www.ricescience.org.

Supplemental Fig. 1. Mean quality distribution of samples in this study.

Supplemental Fig. 2. Length distributions and exon number distributions of lncRNAs and mRNAs.

Supplemental Fig. 3. Expression level of lncRNAs and mRNAs in samples.

Supplemental Table 1. All the primers used for qRT-PCR.

Supplemental Table 2. List of identified 1 655 lncRNA transcripts in DXWR.

Supplemental Table 3. Characteristics of lncRNA transcripts detected in this study.

Supplemental Table 4. Characteristics of mRNAs detected in this study.

Supplemental Table 5. Differentially expressed lncRNAs under drought stress.

Supplemental Table 6. Annotation of the target transcripts of the differentially expressed lncRNAs.

Supplemental Table 7. The target transcripts within 50 kb upstream and downstream of the differentially expressed lncRNAs.

Supplemental Table 8. GO analysis for the target genes of the differentially expressed lncRNAs in the cellular component category.

Supplemental Table 9. GO analysis for the target genes of the differentially expressed lncRNAs in the molecular function category.

Supplemental Table 10. GO analysis for the target genes of the differentially expressed lncRNAs in the biological process category.

Supplemental Table 11. KEGG pathways for the target genes of the differentially expressed lncRNAs.

Ahmad I, Devonshire J, Mohamed R, Schultze M, Maathuis F J M. 2016. Overexpression of the potassium channel TPKb in small vacuoles confers osmotic and drought tolerance to rice., 209(3): 1040–1048.

Anupama A, Bhugra S, Lall B, Chaudhury S, Chugh A. 2018. Assessing the correlation of genotypic and phenotypic responses ofrice varieties under drought stress., 127: 343–354.

Cagirici H B, Alptekin B, Budak H. 2017. RNA sequencing and co-expressed long non-coding RNA in modern and wild wheats., 7(1): 10670.

Cao J N. 2014. The functional role of long non-coding RNAs and epigenetics., 16: 11.

Dr?ge-Laser W, Snoek B L, Snel B, Weiste C. 2018. ThebZIP transcription factor family: An update., 45: 36–49.

Duan H, Tong H, Liu Y Q, Xu Q F, Ma J, Wang C M. 2019. Research advances in the effect of heat and drought on rice and its mechanism., 33(3): 206–218. (in Chinese with English abstract)

Heo J B, Lee Y S, Sung S. 2013. Epigenetic regulation by long noncoding RNAs in plants., 21: 685–693.

Huang L Y, Zhang F, Zhang F, Wang W S, Zhou Y L, Fu B Y, Li Z K. 2014. Comparative transcriptome sequencing of tolerant rice introgression line and its parents in response to drought stress., 15(1): 1026.

Jain M, Nijhawan A, Tyagi A K, Khurana J P. 2006. Validation of housekeeping genes as internal control for studying gene expression in rice by quantitative real-time PCR., 345(2): 646–651.

Jisha V, Dampanaboina L, Vadassery J, Mith?fer A, Kappara S, Ramanan R. 2015. Overexpression of an AP2/ERF type transcription factorconfers biotic and abiotic stress tolerance in rice., 10(6): e0127831.

Joshi R, Wani S H, Singh B, Bohra A, Dar Z A, Lone A A, Pareek A, Singla-Pareek S L. 2016. Transcription factors and plants response to drought stress: Current understanding and future directions., 7: 1029.

Kim D, Langmead B, Salzberg S L. 2015. HISAT: A fast spliced aligner with low memory requirements., 12(4): 357–360.

Kiranmai K, Lokanadha Rao G, Pandurangaiah M, Nareshkumar A, Amaranatha Reddy V, Lokesh U, Venkatesh B, Anthony Johnson A M, Sudhakar C. 2018. A novel WRKY transcription factor,(Lam. Verdc.) enhances drought stress tolerance in transgenic groundnut (L.) plants., 9: 346.

Kong L, Zhang Y, Ye Z Q, Liu X Q, Zhao S Q, Wei L P, Gao G. 2007. Cpc: Assess the protein-coding potential of transcripts using sequence features and support vector machine., 35: W345–W349.

Kroll K W, Mokaram N E, Pelletier A R, Frankhouser D E, Westphal M S, Stump P A, Stump C L, Bundschuh R, Blachly J S, Yan P. 2014. Quality control for RNA-Seq (QuaCRS): An integrated quality control pipeline., 13: 7–14.

Lee D K, Kim H I, Jang G, Chung P J, Jeong J S, Kim Y S, Bang S W, Jung H, Choi Y D, Kim J K. 2015. The NF-YA transcription factorconfers drought stress tolerance of rice in an abscisic acid independent manner., 241: 199–210.

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 S, Scanlon M J, Yu J M, Schnable P S, Timmermans M C P, Springer N M, Muehlbauer G J. 2014. Genome-wide discovery and characterization of maize long non-coding RNAs., 15: R40.

Li S J, Castillo-González C, Yu B, Zhang X R. 2017. The functions of plant small RNAs in development and in stress responses., 90: 654–670.

Li S X, Yu X, Lei N, Cheng Z H, Zhao P J, He Y K, Wang W Q, Peng M. 2017. Corrigendum: Genome-wide identification and functional prediction of cold and/or drought-responsive lncRNAs in cassava., 7: 46795.

Liang Y S, Zheng J, Yan C, Li X X, Liu S F, Zhou J J, Qin X J, Nan W B, Yang Y Q, Zhang H M. 2018. Locating QTLs controlling overwintering trait in Chinese perennial Dongxiang wild rice., 293(1): 81–93.

Liu H P, Yu H Y, Tang G L, Huang T B. 2018. Small but powerful: Function of microRNAs in plant development., 37(3): 515–528.

Lu X K, Chen X G, Mu M, Wang J J, Wang X G, Wang D L, Yin Z J, Fan W L, Wang S, Guo L X, Ye W W. 2016. Genome-wide analysis of long noncoding RNAs and their responses to drought stress in cotton (L.)., 11(6): e0156723.

Luo L J. 2010. Breeding for water-saving and drought-resistance rice (WDR) in China., 61(13): 3509–3517.

Luo X, Bai X, Zhu D, Li Y, Ji W, Cai H, Wu J, Liu B H, Zhu Y M. 2012. GsZFP1, a new Cys2/His2-type zinc-finger protein, is a positive regulator of plant tolerance to cold and drought stress., 235(6): 1141–1155.

Ma J H, Gao X L, Liu Q, Shao Y, Zhang D J, Jiang L N, Li C X. 2017. Overexpression ofincreases drought tolerance through inducing stomatal closure in., 8: 2036.

Mao X Z, Cai T, Olyarchuk J G, Wei L P. 2005. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary., 21(19): 3787–3793.

Menguer P K, Sperotto R A, Ricachenevsky F K. 2017. A walk on the wild side:species as source for rice abiotic stress tolerance., 40: 238–252.

Moran E, Lauder J, Musser C, Stathos A, Shu M J. 2017. The genetics of drought tolerance in conifers., 216(4): 1034–1048.

Munné-Bosch S, Falara V, Pateraki I, López-Carbonell M, Cela J, Kanellis A K. 2009. Physiological and molecular responses of the isoprenoid biosynthetic pathway in a drought-resistant Mediterranean shrub,exposed to water deficit., 166(2): 136–145.

Nawaz F, Naeem M, Zulfiqar B, Akram A, Ashraf M Y, Raheel M, Shabbir R N, Hussain R A, Anwar I, Aurangzaib M. 2017. Understanding brassinosteroid-regulated mechanisms to improve stress tolerance in plants: A critical review., 24(19): 15959–15975.

Pauli A, Valen E, Lin M F, Garber M, Vastenhouw N L, Levin J Z, Fan L, Sandelin A, Rinn J L, Regev A, Schier A F. 2012. Systematic identifcation of long noncoding RNAs expressed during zebrafsh embryogenesis., 22(3): 577–591.

Pertea M, Pertea G M, Antonescu C M, Chang T C, Mendell J T, Salzberg S L. 2015. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads., 33(3): 290–295.

Phukan U J, Jeena G S, Shukla R K. 2016. WRKY transcription factors: Molecular regulation and stress responses in plants., 7: 760.

Punta M, Coggill P C, Eberhardt R Y, Mistry J, Tate J, Boursnell C, Pang N, Forslund K, Ceric G, Clements J, Heger A, Holm L, Sonnhammer E L L, Eddy S R, Bateman A, Finn R D. 2012. The Pfam protein families database., 40: D290– D301.

Qin T, Zhao H Y, Cui P, Albesher N, Xiong L M. 2017. A nucleus-localized long non-coding RNA enhances drought and salt stress tolerance., 175(3): 1321–1336.

Seo J S, Joo J, Kim M J, Kim Y K, Nahm B H, Song S I, Cheong J J, Lee J S, Kim J K, Choi Y D. 2011. OsbHLH148, a basic helix-loop-helix protein, interacts with OsJAZ proteins in a jasmonate signaling pathway leading to drought tolerance in rice., 65(6): 907–921.

Shanker A K, Maheswari M, Yadav S K, Desai S, Bhanu D, Attal N B, Venkateswarlu B. 2014. Drought stress responses in crops., 14(1): 11–22.

Song X W, Cao X F. 2017. Transposon-mediated epigenetic regulation contributes to phenotypic diversity and environmental adaptation in rice., 36: 111–118.

Song Y, Jing S J, Yu D Q. 2009. Overexpression of the stress- inducedimproves osmotic stress tolerance in., 54(24): 4671–4678.

Song Y, Chen L G, Zhang L P, Yu D Q. 2010. Overexpression ofgene interferes in the abscisic acid signal and auxin transport pathway of., 35(3): 459–471.

Sun L, Luo H T, Bu D C, Zhao G G, Yu K T, Zhang C H, Liu Y N, Chen R S, Zhao Y. 2013. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts., 41: e166.

Sun X, Zhang Y P, Zhu X D, Korir N K, Tao R, Wang C, Fang J G. 2014. Advances in identification and validation of plant microRNAs and their target genes., 152(2): 203–218.

Trapnell C, Williams B A, Pertea G, Mortazavi A, Kwan G, van Baren M J, Salzberg S L, Wold B J, Pachter L. 2010. Transcript assembly and quanti?cation by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation., 28: 511–515.

Wang L G, Park H J, Dasari S, Wang S Q, Kocher J P, Li W. 2013. Cpat: Coding-potential assessment tool using an alignment-free logistic regression model., 41(6): e74.

Wang L K, Feng Z X, Wang X, Wang X W, Zhang X G. 2010. Degseq: An R package for identifying differentially expressed genes from RNA-seq data., 26(1): 136–138.

Wang L L, Yu C C, Chen C, He C L, Zhu Y G, Huang W C. 2014. Identification of rice Di19 family revealsinvolved in drought resistance., 33(12): 2047–2062.

Xiong L Z, Yang Y N. 2003. Disease resistance and abiotic stress tolerance in rice are inversely modulated by an abscisic acid- inducible mitogen-activated protein kinase., 15(3): 745–759.

Xu D Q, Huang J, Guo S Q, Yang X, Bao Y M, Tang H J, Zhang H S. 2008. Overexpression of a TFIIIA-type zinc finger protein geneenhances drought and salt tolerance in rice (L.)., 582(7): 1037–1043.

Yang J, Chen X R, Zhu C L, Peng X S, He X P, Fu J R, Ouyang L J, Bian J M, Hu L F, Sun X T, Xu J, He H H. 2015. RNA-seq reveals differentially expressed genes of rice () spikelet in response to temperature interacting with nitrogen at meiosis stage., 16: 959.

Yasuda K, Ito M, Sugita T, Tsukiyama T, Saito H, Naito K, Teraishi M, Tanisaka T, Okumoto Y. 2013. Utilization of transposable element as a novel genetic tool for modification of the stress response in rice., 32(3): 505–516.

Ye N H, Zhu G H, Liu Y G, Li Y X, Zhang J H. 2011. ABA controls H2O2accumulation through the induction ofin rice leaves under water stress., 52(4): 689–698.

Yin M Z, Wang Y P, Zhang L H, Li J Z, Quan W L, Yang L, Wang Q F, Chan Z L. 2017. TheCys2/His2 zinc finger transcription factor ZAT18 is a positive regulator of plant tolerance to drought stress., 68(11): 2991–3005.

Young M D, Wakefield M J, Smyth G K, Oshlack A. 2010. Gene ontology analysis for RNA-seq: Accounting for selection bias., 11(2): R14.

Yu X, Niu X L, Yang S H, Li Y X, Liu L L, Tang W, Liu Y S. 2011. Research on heat and drought tolerance in rice (L.) by overexpressing transcription factor., 44(20): 4142–4149.

Yuan J P, Li J R, Yang Y, Tan C, Zhu Y M, Hu L, Qi Y J, Lu Z J. 2018. Stress-responsive regulation of long non-coding RNA polyadenylation in., 93(5): 814–827.

Zaynab M, Fatima M, Abbas S, Umair M, Sharif Y, Raza M A. 2018. Long non-coding RNAs as molecular players in plant defense against pathogens., 121: 277–282.

Zhang F T, Luo X D, Zhou Y, Xie J K. 2016. Genome-wide identification of conserved microRNA and their response to drought stress in Dongxiang wild rice (Griff.)., 38(4): 711–721.

Zhang F T, Luo Y, Zhang M, Zhou Y, Chen H P, Hu B L, Xie J K. 2018. Identification and characterization of drought stress- responsive novel microRNAs in Dongxiang wild rice., 25(4): 175–184.

Zhang G Y, Duan A G, Zhang J G, He C Y. 2017. Genome-wide analysis of long non-coding RNAs at the mature stage of sea buckthorn (Linn) fruit., 596: 130–136.

Zhang L Y, Gu L K, Ringler P, Smith S, Rushton P J, Shen Q J. 2015. Three WRKY transcription factors additively repress abscisic acid and gibberellin signaling in aleurone cells., 236: 214–222.

Zhang W, Han Z X, Guo Q L, Liu Y, Zheng Y X, Wu F L, Jin W B. 2014. Identification of maize long non-coding RNAs responsive to drought stress., 9(6): e98958.

Zhang X, Zhou S X, Fu Y C, Su Z, Wang X K, Sun C Q. 2006. Identification of a drought tolerant introgression line derived from Dongxiang common wild rice (.Griff.)., 62: 247–259.

Zhao J F, Favero D S, Qiu J W, Roalson E H, Neff M M. 2014. Insights into the evolution and diversification of the AT-hook motif nuclear localized gene family in land plants., 14: 266.

Zhou L G, Liu Z C, Liu Y H, Kong D Y, Li T F, Yu S W, Mei H W, Xu X Y, Liu H Y, Chen L, Luo L J. 2016. A novel geneimproves both drought avoidance and drought tolerance in rice., 6: 30264.

Zhu B Z, Yang Y F, Li R, Fu D Q, Wen L W, Luo Y B, Zhu H L. 2015. RNA sequencing and functional analysis implicate the regulatory role of long non-coding RNAs in tomato fruit ripening., 66(15): 4483–4495.

Copyright ? 2020, China National Rice Research Institute. Hosting by Elsevier B V

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/)

Peer review under responsibility of China National Rice Research Institute

http://dx.doi.org/10.1016/j.rsci.2019.12.003

20 July 2018;

27 September 2018

s:Zhang Fantao (zhang84004@163.com); Xie Jiankun (xiejiankun11@163.com)

(Managing Editor: Li Guan)