Genome-wide analysis of long non-coding RNAs (lncRNAs) in two contrasting rapeseed (Brassica napus L.) genotypes subjected to drought stress and re-watering

Background Drought stress is a major abiotic factor that affects rapeseed (Brassica napus L.) productivity. Though previous studies indicated that long non-coding RNAs (lncRNAs) play a key role in response to drought stress, a scheme for genome-wide identification and characterization of lncRNAs’ response to drought stress is still lacking, especially in the case of B. napus. In order to further understand the molecular mechanism of the response of B. napus to drought stress, we compared changes in the transcriptome between Q2 (a drought-tolerant genotype) and Qinyou8 (a drought-sensitive genotype) responding drought stress and rehydration treatment at the seedling stage. Results A total of 5546 down-regulated and 6997 up-regulated mRNAs were detected in Q2 compared with 7824 and 10,251 in Qinyou8, respectively; 369 down-regulated and 108 up- regulated lncRNAs were detected in Q2 compared with 449 and 257 in Qinyou8, respectively. LncRNA-mRNA interaction network analysis indicated that the co-expression network of Q2 was composed of 145 network nodes and 5175 connections, while the co-expression network of Qinyou8 was composed of 305 network nodes and 22,327 connections. We further identified 34 transcription factors (TFs) corresponding to 126 differentially expressed lncRNAs in Q2, and 45 TFs corresponding to 359 differentially expressed lncRNAs in Qinyou8. Differential expression analysis of lncRNAs indicated that up- and down-regulated mRNAs co-expressed with lncRNAs participated in different metabolic pathways and were involved in different regulatory mechanisms in the two genotypes. Notably, some lncRNAs were co-expressed with BnaC07g44670D, which are associated with plant hormone signal transduction. Additionally, some mRNAs co-located with XLOC_052298, XLOC_094954 and XLOC_012868 were mainly categorized as signal transport and defense/stress response. Conclusions The results of this study increased our understanding of expression characterization of rapeseed lncRNAs in response to drought stress and re-watering, which would be useful to provide a reference for the further study of the function and action mechanisms of lncRNAs under drought stress and re-watering.


Background
Drought is one of the vital factors limiting crop productivity and survival. Due to the ongoing global climate change, more and more research has focused on understanding the mechanisms of how crops resist drought stress and improve their resistance level [1][2][3][4][5][6][7][8]. Plants sense drought signals and produce second messenger substances, such as Ca 2+ , phosphatidylinositol and reactive oxygen species (ROS) [9,10], while causing an increase in intracellular calcium ion concentration, initiating a cascade network of protein phosphorylation pathways. Finally, the target proteins are directly involved in the protection of cells, or regulate the expression of a series of specific stress-related genes through TFs (MYC/MYB, ABF, CBF/DREB, bZIP, etc.), thereby protecting the cells and improving the resistance of plants to adversity [11][12][13]. Although rapid developments in modern molecular biology have gradually uncovered the molecular mechanisms of plant drought resistance, developing drought-resistant plants to cope with drought stress will remain a substantive challenge in the future.
Long non-coding RNA (lncRNA) is a type of RNA transcripts which is more than 200 nucleotides in length and has no or limited protein coding abilities [14][15][16]. A growing body of evidence has shown that lncRNAs exert their regulatory effects on gene expression levels, involving epigenetic regulation, transcriptional regulation, and posttranscriptional regulation in the form of RNA [17][18][19][20][21][22][23][24][25]. With the advantage of next-generation sequencing technologies and bioinformatics approaches, many lncRNAs have been discovered in model plants, such as Arabidopsis [26][27][28][29], wheat [30], maize [31][32][33] and rice [34], indicating that lncRNAs play an important role in various biological processes of plant development and stress response. Recent research has confirmed that lncRNAs respond to abiotic stresses [31,35,36], including drought stress. For example, 664 drought-responsive lncRNAs were analyzed in maize [31]. Under drought stress, 2542 lncRNA candidates have been identified from Populus trichocarpa, 504 of which were found to be drought-responsive [37]. In Arabidopsis, 1832 lncRNAs changed after 2 h and/or 10 h of drought, cold, high-salt, and/or abscisic acid (ABA) treatments [29]. In maize, 664 transcripts were confirmed as drought-responsive lncRNAs, 8 out of which were proved as precursors of miRNAs [31]. In rice, pre-miRNA expression profiling indicated that miR171f is involved in the progression of rice root development and growth, as well as the response to drought stress [38]. In cotton, long intervening / intergenic noncoding RNAs (lincRNAs) XLOC 063105 and XLOC 115463, were involved in drought stress response by regulating neighboring genes [39]. Furthermore, 19 lncRNAs (17 lincRNAs and 2 natural antisense transcripts (NATs)) in foxtail millet responded to polyethylene glycol-6000 (PEG)-induced drought stress, only one of the drought-responsive lncRNA had synteny with its sorghum counterpart [40]. Qin et al. (2017) identified an Arabidopsis lncRNA, drought-induced lncRNA (DRIR), which responds to drought and salt stress. DRIR can be significantly activated by drought and salt stress as well as by abscisic acid (ABA) treatment [41]. In addition, in cassava, 318 lncRNAs were identified, which were responsive to cold and/or drought stress, and which are associated with hormone signal transduction, biosynthesis of secondary metabolites, and the sucrose metabolism pathway [42]. Additionally, numerous lncRNAs involved in the regulation of gene expression in response to stress have been identified and characterized in Brassica [43][44][45][46]. In Chinese cabbage (Brassica rapa ssp. chinensis), 4594 putative lncRNAs were identified to response to heat stress, 25 of which were co-expressed with 10 heat responsive genes [47]. In Brassica rapa L., 549 lncRNAs were identified significantly altered their expression in response to cold treatment, and short-term cold treatment induced natural antisense transcripts (NATs) in BrFLC and BrMAF genes which are involved in vernalization were identified [48]. Summanwar et al. (2019) identified 530 differentially expressed lncRNAs from the roots of clubroot-susceptible and -resistant Brassica napus lines. Twenty-four differentially expressed lncRNAs were identified from chromosome A08 which has been reported to confer resistance to different P. brassicae pathotypes [49]. In Brassica juncea, 1614 differentially expressed lncRNAs response to heat and drought stress, and some lncRNAs were co-expressed with TFs which are involved in abiotic stress response [50].
Rapeseed (Brassica napus L.) is an important oilseed crop in the world [51]. It is vulnerable to drought, which influences the production of rapeseed substantially [52][53][54]. Although many lncRNAs have been found in different plant species, indicating that lncRNAs can play an important role in response to abiotic stresses, a genome-wide identification and characterization of responses of lncRNAs to drought stress and rehydration treatments is still lacking, especially in B. napus. In order to further understand the molecular mechanisms of the response of B. napus to drought stress and re-watering, we compared changes in transcriptome between Q2 (a drought-tolerant genotype) and Qinyou8 (a drought-sensitive genotype) in response to drought stress and rehydration treatments at the seedling stage, and identified the lncRNAs involved in drought stress and rehydration treatments. The present study used a co-expression-based method, in which lncRNA functions were predicted, based on the functions of their co-expressed protein-coding genes [55]. Therefore, the lncRNA-mRNA co-expression network was constructed for pathway enrichment analysis. Moreover, the lncRNA-mRNA co-expression network of plant hormone signal transduction was analyzed to further explore the potential roles of differentially expressed lncRNAs in response to drought stress and re-watering.

Results
Phenotypes of rapeseed seedlings under drought stress (DS) and re-watering (RW) treatments Rapeseed seedlings responded differently under the DS and RW treatments (Fig. 1). The fresh weight of Q2 under DS reached 70.43% of well-watered (WW), which was significantly higher than that of Qinyou8. Additionally, the fresh weight of Q2 under RW reached 82.76% of WW, which was significantly higher than Qinyou8. Therefore, we can see that the recovery ability of Q2 is better than that of Qinyou8 after re-watering.

Differentially expressed lncRNAs and mRNAs under drought stress and re-watering
In this study, RNAs were extracted from 12 samples (two treatments, two test materials, three biological replicates) and tested their quality before performing RNA sequencing (Additional file 1). We acquired clean reads by removing low-quality reads from the RNA-seq data. The QC and GC contents were calculated from clean data to assess the quality of the sequencing data (Additional file 2). The clean datasets were mapped to the Brassica napus L. genome. All results indicated that the RNA-seq data were very reliable. The expression level of all transcripts, including lncRNAs and mRNAs, were identified using FPKM, which was systematically estimated and the differential transcript analysis done using cuffdiff with a threshold of q value < 0.05. Compared with the expression of lncRNAs in drought stress, 477 lncRNAs (369 down-regulated, 108 up-regulated) of Q2 and 706 lncRNAs (449 down-regulated, 257 upregulated) of Qinyou8 were differentially expressed after re-watering (Fig. 2). In addition, there were 12,543 mRNAs (5546 down-regulated, 6997 up-regulated) and 18,075 mRNAs (7824 down-regulated, 10,251 upregulated) differentially expressed in Q2 and Qinyou8, respectively (Fig. 2).

qRT-PCR validation
To validate the expression data from RNA-seq, nine lncRNAs which were differentially expressed in two genotypes were selected for real-time RT-PCR analysis. The squared of the pearson's correlation coefficient (the determinant coefficient) of lncRNAs expression level was calculated. As shown in Fig. 3, the lncRNA expression level using RNA-seq was significantly (R 2 = 0.91519, slope = 0.91646) correlated with those using qRT-PCR. For example, the relative expression of XLOC_012868 was increased in Q2 but decreased in Qinyou8, which was consistent with the RNA-seq result (Additional file 3). The real-time PCR results verify the expression patterns obtained with transcriptome sequencing, indicating that the lncRNAs expression profile based on RNA-seq data is reliable.

Functions of differentially expressed lncRNAs based on lncRNA-mRNA co-expression network
To further characterize the role of differentially expressed lncRNAs, we used the lncRNA-mRNA relationship pairs to construct interactive networks. Coexpression network analysis indicated that the coexpression network of Q2 was composed of 145 network nodes and 5175 connections, while the co-expression network of Qinyou8 was composed of 305 network Comparisons of fresh weight among the treatments. Experiments were repeated three times and vertical bars indicate standard errors. "*" indicates the significance of the difference at the 0.05 level. WW = well-watered; DS = drought stress; RW = re-watering nodes and 22,327 connections. In Q2, there were 5175 lncRNA-mRNA pairs, including 1481 mRNAs and 145 lncRNAs, respectively (Additional file 4). Similarly, there were 22,327 lncRNA-mRNA pairs in Qinyou8, which included 3200 mRNAs and 305 lncRNAs, respectively (Additional file 4). The lncRNA-mRNA pairs with the same expression trend were much more than those with the opposite expression trend in two genotypes. There were 6 and 4 opposite-trend pairs in Q2 and Qinyou8, respectively, which suggested the candidate lncRNAs function in drought and re-watering processes (Additional file 5). It has been shown that one lncRNA may regulate multiple protein-coding genes, and vice versa [56][57][58]. From the co-expression network of Q2, we know that 1 mRNA may correlate with 1 to 18 lncRNAs, and that 1 lncRNA may correlate with 1 to 375 mRNAs. Moreover, the co-expression network of Qinyou8 indicated that 1 mRNA may correlate with 1 to 43 lncRNAs, and 1 lncRNA may correlate with 1 to 375 mRNAs. XLOC_071559 was the largest node in the network in both genotypes, respectively.
Studies have shown that lncRNA can indirectly affect the expression of mRNA, and can also directly bind to mRNA, thus affecting translation [59][60][61], shearing [62,63], and degradation of mRNA [64]. Currently, the mechanism of interaction between lncRNA and mRNA has not yet become clear. To reveal potential functions of the differentially expressed lncRNAs under drought stress and re-watering, we analyzed Gene Ontology (GO) terms of target genes of differentially expressed lncRNAs. This analysis was performed to determine the major molecular functions, biological processes, and cellular components with which the target genes of differentially expressed lncRNAs were associated.
The down-regulated mRNAs, co-expressed with differentially expressed lncRNAs, were assigned to 32 and 34 significant terms in Q2 and Qinyou8, respectively (Fig. 4a). For down-regulated mRNAs co-expressed with differentially expressed lncRNAs in Q2, the most significant GO terms for biological process were oxidationreduction process (GO:0055114), protein dephosphorylation (GO: 0006470), dephosphorylation (GO:0016311), Fig. 2 The numbers of differentially expressed lncRNAs and mRNAs in two genotypes (Q2 and Qinyou 8) in response to drought stress and re-watering treatments response to abiotic stimulus (GO:0009628), response to water stimulus (GO:0009415) and sucrose metabolic process (GO:0005985). As far as molecular functions are concerned, nucleic acid binding transcription factor activity (GO:0001071), sequence-specific DNA binding transcription factor activity (GO:0003700), cofactor binding (GO:0048037), sequence-specific DNA binding (GO:0043565), phosphoric ester hydrolase activity (GO: 0042578), protein serine/threonine phosphatase activity (GO:0004722) and phosphoprotein phosphatase activity (GO:0004721) were the important significantly enriched GO terms. The GO terms of transcription factor complex (GO:0005667) and CCAAT-binding factor complex (GO:0016602) were the most important significant terms for cellular components. In Qinyou8, for down-regulated mRNAs co-expressed with differentially expressed lncRNAs, the important GO terms for biological process were single-organism metabolic process (GO:0044710), oxidation-reduction process (GO:0055114), protein dephosphorylation (GO:0006470), response to abiotic stimulus (GO:0009628), response to water stimulus (GO: 0009415) and protein serine/threonine phosphatase activity (GO:0004722). As far as molecular functions are concerned, three GO terms, namely, oxidoreductase activity (GO: 0016491), nucleic acid binding transcription factor activity (GO:0001071) and sequence-specific DNA binding transcription factor activity (GO:0003700), demonstrated significant enrichment. With respect to cellular components, transcription factor complex (GO: 0005667) and CCAAT-binding factor complex (GO: 0016602) were the most significantly enriched GO terms.
The up-regulated mRNAs co-expressed with differentially expressed lncRNAs were assigned to 23 and 31 significant GO terms in Q2 and Qinyou8, respectively (Fig. 4b). Of the enriched GO terms in the biological process category for up-regulated mRNAs co-expressed with differentially expressed lncRNAs in Q2, primary metabolic process (GO:0044238) was the most dominant group, followed by organic substance metabolic process (GO:0071704), macromolecule metabolic process (GO: 0043170), protein metabolic process (GO:0019538), cellular macromolecule metabolic process (GO: 0044260) and cellular protein metabolic process (GO:0044267). Among the molecular functions, structural molecule activity (GO:0005198) and structural constituent of ribosome (GO:0003735) were the most dominant groups in Q2. In the cellular component category, the significant terms were cytoplasm (GO:0005737), cytoplasmic part (GO:0044444), ribosome (GO:0005840), ribonucleoprotein complex (GO:0030529) and translation (GO: 0006412). Additionally, in Qinyou8, the GO terms of , ribosome (GO:0005840) and photosystem II (GO: 0009523) were the dominant groups. These findings suggest that stress-responsive lncRNAs may regulate genes involved in many biological processes, including signal transduction, energy synthesis, molecule metabolism, transcription and translation, in response to drought stress and re-watering.
We also analyzed the statistical enrichment of the mRNAs co-expressed with differentially expressed lncRNAs in KEGG. There were 18 and 18 KEGG Fig. 4 Gene Ontology (GO) classifications of the co-expressed mRNAs of the differentially expressed lncRNAs. The mRNAs co-expressed with lncRNAs are divided into three main categories by GO analysis: biological process, molecular function and cellular component. The x-axis indicates the number of genes in a sub-category, and the y-axis indicates the sub-categories. a The most highly enriched GO terms for the down-regulated genes in Q2 and Qinyou8. b The most highly enriched GO terms for the up-regulated genes in Q2 and Qinyou8 pathways identified significantly in Q2 and Qinyou8, respectively, using pathway enrichment analysis (p < 0.05). KEGG analysis showed that there were 19 pathways identified that significantly related to downregulated mRNAs co-expressed with differentially expressed lncRNAs of Q2 (Fig. 5a), including plant hormone signal transduction (ko04075), glycolysis/gluconeogenesis (ko00010), fatty acid metabolism (ko01212), valine, leucine and isoleucine degradation (ko00280), alanine, aspartate and glutamate metabolism (ko00250) and arginine and proline metabolism (ko00330). Moreover, 3 identified pathways significantly related to upregulated mRNAs co-expressed with differentially expressed lncRNAs of Q2 (Fig. 5b), including ribosome (ko03008), carbon fixation in photosynthetic organisms (ko00710), and pyruvate metabolism (ko00620). In Qinyou8, 17 pathways were identified that were significantly related to down-regulated mRNAs co-expressed with differentially expressed lncRNAs and 7 identified pathways significantly related to up-regulated mRNAs Fig. 5 KEGG pathways analysis. Top 20 pathways for the co-expressed mRNAs of the differentially expressed lncRNAs. The y-axis corresponds to the KEGG pathway with a q value ≤0.05, and the x-axis shows the enrichment ratio between the number of differentially expressed genes and all UniGenes enriched in a particular pathway. The color of the dot represents q value, and the size of the dot represents the number of differentially expressed genes mapped to the reference pathways. a KEGG pathway classification of the down-regulated mRNAs co-expressed with lncRNAs in Q2. b KEGG pathway classification of the up-regulated mRNAs co-expressed with lncRNAs in Q2. c KEGG pathway classification of the down-regulated mRNAs co-expressed with lncRNAs in Qinyou8. d KEGG pathway classification of the up-regulated mRNAs co-expressed with lncRNAs in Qinyou8 co-expressed with differentially expressed lncRNAs. The most down-regulated mRNAs co-expressed with lncRNAs of Qinyou8 were significantly enriched for protein processing in endoplasmic reticulum (ko04141), fatty acid metabolism (ko01212), fatty acid degradation (ko00071), alanine, aspartate and glutamate metabolism (ko00250), galactose metabolism (ko00052), as well as arginine and proline metabolism (ko00330) (Fig. 5c), while the most up-regulated mRNAs co-expressed with lncRNAs of Qinyou8 denoted their involvement in ribosome (ko03010), photosynthesis (ko00195), and photosynthesis -antenna proteins (ko00196) (Fig. 5d).

Discussion
Several recent studies have revealed that lncRNAs play an important role in response to drought stress [31,39,41,66]. Accordingly, we constructed lncRNA and mRNA libraries, and annotated, identified, and verified those lncRNAs that are involved in drought stress and re-watering.
Differential mRNAs and lncRNAs expression in two contrasting genotypes under drought stress and rewatering In our research, we systematically identified and analyzed B. napus mRNAs and lncRNAs, which respond to drought stress and rehydration. In the comparison groups with two different genotypes, 5546 downregulated and 6997 up-regulated mRNAs were detected in Q2 compared to 7824 and 10,251 in Qinyou8, respectively; 369 down-regulated and 108 up-regulated lncRNAs were detected in Q2 compared to 449 and 257 in Qinyou8, respectively. Interestingly, we found that there were 229 lncRNAs (169 down-regulated, 44 upregulated) in both genotypes, among which, 1 lncRNA XLOC_012868 was up-regulated in the drought-tolerant genotype and down-regulated in drought-susceptible genotype; conversely, 15 lncRNAs were down-regulated in the drought-tolerant genotype and up-regulated in the drought-susceptible genotype (Fig. 7). From the above, we know that the response of these two genotypes is different under drought stress and rehydration conditions. In Qinyou8, the number of differentially expressed mRNAs and lncRNAs was significantly higher than Q2.
Altered splicing is one of the mechanisms for lncRNA transcripts to affect gene expression in many physiological processes [67][68][69]. In Q2, 477 lncRNA transcripts from 469 lncRNA genes were identified, in which 8 lncRNA coding genes were alternatively spliced. Similarly, in Qinyou8, 706 lncRNA transcripts from 688 lncRNA genes were identified, in which 18 lncRNA coding genes were alternatively spliced (Additional file 7). These alternately spliced lncRNA coding genes may be involved in drought and re-watering processes. Additionally, 9 identified lncRNAs were chosen for qRT-PCR validation, and the results confirmed the sequencing results.

Differentially expressed lncRNAs specifically enriched in GO and KEGG pathways
With advances in next-generation sequencing technology, many investigations have shown that lncRNAs exert their regulatory effects on gene expression levels, including epigenetic regulation, transcriptional regulation, and post-transcriptional regulation in the form of RNA [19]. It is known that sequence-specific DNA binding transcription factor activity [42,70], response to stimulus [71], response to abiotic stimulus [70], biosynthetic process [70], structural constituent of ribosome [58], photosynthesis [72] and oxidoreductase activity [72], which are regulated by some lncRNAs, have been reported in response to abiotic stresses, and these GO terms were identified in this study. To determine the similarity and differences between the two genotypes, the significantly enriched GO terms were compared. In our study, there were more significant GO terms in Qinyou8 than Q2 under drought stress and re-watering, indicating that there were differences in responses to drought stress and re-watering between the two genotypes. We found that phosphoprotein phosphatase activity, protein metabolic process, and sequence-specific DNA binding were significantly and specially enriched in Q2, while single-organism metabolic process, photosynthesis, and oxidoreductase activity were significantly and specially enriched in Qinyou8. Additionally, lncRNAs have been recognized as powerful regulators of pathways in response to drought stress, including ribosome, photosynthesis [73], and plant hormone signal transduction [42,74]. The ribosome pathway was simultaneously significant in both genotypes, and the differential lncRNA target genes were up-regulated in this pathway. It is worth noting that the pathway of plant hormone signal transduction was significantly and specially enriched in Q2, a total of 36 mRNAs co-expressed with 41 lncRNAs were assigned to plant hormone signal transduction. Furthermore, many down-regulated mRNAs co-expressed with lncRNAs involved in protein processing in the endoplasmic reticulum, and upregulated mRNAs co-expressed with lncRNAs belonging to photosynthesis were significantly and specially enriched in Qinyou8. A total of 7 and 5 mRNAs coexpressed with lncRNAs were assigned into photosynthesis and photosynthesis-antenna proteins, respectively. The genes involved in photosynthesis were generally down-regulated by drought [75,76]. Compared with the DS treatment, photosynthesis (ko00195) and photosynthetic antenna protein (ko00196) pathways were significantly enriched in Qinyou8 under the RW treatment, indicating that the short-term drought stress did not cause significantly damage to photosynthesis of Q2, but did some damage to Qinyou8. In Qinyou8, the genes involved in photosynthesis (ko00195) and photosynthetic antenna protein (ko00196) were up-regulated to restore normal photosynthesis and thus restore growth. These results indicate that lncRNAs could play a role in many biological processes responding to drought stress and rewatering through regulating gene network, and that up-and down-regulated mRNAs co-expressed with lncRNAs participate in different metabolic pathways and are involved in different regulation mechanisms. Taken together, our results suggest that the two different genotypes implement divergent mechanisms to modulate the response to drought stress and re-watering treatment.
Analysis of plant signal transduction using lncRNA-mRNA co-expression network Regulation on the co-expression network may be the possible mechanisms in response to stress for lncRNAs [18,31]. Although a large number of lncRNAs were identified to be related with many biological processes, a limited number of lncRNAs were screened out to contribute to plant hormone signal transduction by using lncRNA-mRNA co-expression analysis. In Q2, the coexpression network of plant hormone signal transduction contained 157 matched lncRNA-mRNA pairs, including 41 lncRNAs and 36 mRNAs ( Fig. 8a and Additional file 8). The co-expression network of plant hormone signal transduction of Qinyou8 was composed of 120 lncRNAs and 51 mRNAs with 352 matched lncRNA-mRNA pairs ( Fig. 8b and Additional file 8). The lncRNAs involved in plant hormone signal transduction had the same expression direction with the target genes in two genotypes, proving the expression of lncRNAs promoted the function of the target genes. In this pathway, target genes of differentially expressed lncRNAs were involved in auxin, cytokinin, gibberellin, and abscisic acid signaling pathways in both genotypes. Some target genes of differentially expressed lncRNAs related to the ethylene and salicylic acid signaling pathways were specifically expressed in Q2, while target genes of differentially expressed lncRNAs involved in the two signaling pathways of BR and jasmonic acid were specifically expressed in Qinyou8. Among these signaling pathways, more of the mRNAs, which co-expressed with differentially expressed lncRNAs, were associated with the ABA signaling pathways than those of other phytohormones, which is consistent with previous studies that had considered ABA to be an early warning signal for plant responses to drought stress [77,78].
Auxin (IAA) as a phytohormone, is essential for signaling, transport, growth and development of a plant [79]. Auxin binds to the TRANSPORT INHIBITOR RESPONSE 1/AUXIN SIGNALLING F-BOX proteins (TIR1/AFBs) and the AUXIN/INDOLE-3-ACETIC ACID (Aux/IAA) proteins. When the level of IAA is low, the Aux/IAA protein forms a heterodimer with the auxin response factor (ARF) to inhibit gene transcription. Conversely, the Aux/IAA protein is degraded, which results in derepression of the ARF transcriptional regulation and expression of the auxin response gene [80]. Currently, IAA early response genes mainly include AUXIN/INDOLE-3-ACETIC ACID (Aux/IAA), Gretchen Hagen 3 (GH3) and Small Auxin-Up RNAs (SAUR), which are auxin-induced primitive expression genes [81]. Among them, Aux/IAA protein plays a very important role in the IAA signal transduction pathway, and it acts as a transcriptional repressor in the signal transduction pathway [82]. The GH3 gene encodes an auxin-binding enzyme that acts as a feedback regulator of auxin by reducing the level of beneficial auxin [83]. In Q2, the co-expression network of IAA signal transduction contained 21 matched lncRNA-mRNA pairs, including 16 lncRNAs and 3 mRNAs. In Qinyou8, the co-expression network of IAA signal transduction contained 56 matched lncRNA-mRNA pairs, which included 46 lncRNAs and 12 mRNAs. Drought stress and re-watering regulated the expression of Aux/IAA (1 differentially expressed mRNA co-expressed with lncRNAs in Q2, and 6 differentially expressed mRNAs co-expressed with lncRNAs in Qinyou8), and GH3 (1 differentially expressed mRNA co-expressed with lncRNAs in Q2, and 3 differentially expressed mRNAs co-expressed with lncRNAs in Qinyou8) genes in the two genotypes. In Q2, down-regulated XLOC_042431, XLOC_071559, XLOC_ 095305, XLOC_100682, XLOC_019521 and XLOC_ 042894, targeting down-regulated BnaC06g05090D (encoding Aux/IAA), possibly take part in regulating the IAA signal transduction pathway in a positive way. Furthermore, down-regulated XLOC_098397, XLOC_034532 and XLOC_038342, targeting down-regulated BnaA05g147 80D (encoding GH3), facilitating the level of beneficial auxin. It is suggested that down regulation of these lncRNAs expression in Q2 led to enhance IAA signal, which may accelerate vegetative growth by cell enlargement. In Qinyou8, up-regulated XLOC_017878, XLOC_ 042549, and XLOC_028678, targeting up-regulated BnaC01g06240D, BnaC03g78000D, and BnaC08g43830D (encoding Aux/IAA), respectively. Additionally, upregulated XLOC_017878, targeting up-regulated BnaA 09g42140D and BnaC08g34560D (encoding GH3). Up regulation of these lncRNAs expression in Qinyou8 led to weakened IAA signal, which may inhibit vegetative growth.
Cytokinin (CK) plays an important role in various physiological functions in plants, such as promoting cell division, inducing shoot formation and promoting its growth [79]. Cytokinin signaling is based on a twocomponent signaling system (TCS), which is mainly composed of Arabidopsis histidine kinases (AHKs), Arabidopsis histidine phosphotransfer proteins (AHPs) and Arabidopsis response regulators (ARRs). Firstly, the cytokinin receptor binds to cytokinin and then to autophosphorylates. Subsequently, it transfers the phosphate group to a phosphotransferase of the cytoplasm through transmembrane transport; the phosphorylated AHPs can then enter the nucleus and transfer the phosphate group to the response regulator, thereby inducing gene expression and regulating plant growth and development [84]. The type-B response regulators (B-ARR) function as positive regulators of cytokinin signaling, while the type-A response regulators (A-ARR) function as a downstream signal that acts as the negative regulators of cytokinin signaling and also inhibits the signal transmission of B-ARR [85]. In Q2, the co-expression network of CK signal transduction contained 7 matched lncRNA-mRNA pairs, including 7 lncRNAs and 3 mRNAs. Down-regulated BnaA01g17750D is involved in encoding B-ARR gene, was targeted by down-regulated XLOC_075476 and XLOC_074677, indicating that down-regulated XLOC_075476 and XLOC_074677 are likely to weakened CK signal, which may inhibit the seedling growth of Q2. In Qinyou8, the co-expression network of IAA signal transduction contained 27 matched lncRNA-mRNA pairs, which included 25 lncRNAs and 5 mRNAs. Down-regulated BnaC06g 18770D is involved in encoding A-ARR gene, was targeted by 9 down-regulated lncRNAs. It is suggested that 9 down-regulated lncRNAs of Qinyou8 are likely to enhance CK signal, which may benefit rapeseed seedling growth.
Gibberellin (GA) plays an important role in all stages of plant growth and development, and it participates in various physiological processes that regulate plant growth and development. One of the most significant effects is the promotion of internode elongation, which promotes plant growth [86]. GIBBERELLIN INSENSI-TIVE DWARF1 (GID1) receptor is a soluble protein that is localized to both cytoplasm and nucleus. GID1 protein can specifically bind to active GA and further bind with DELLA protein to form GID1-GA-DELLA [87]. By mediating the degradation of or inhibiting the activity of DELLA protein, the GID1-GA-DELLA disinhibits DELLA protein from the GA reactive system, and then activates the GA reactive gene [88]. When GA is at a low level, GID1 does not bind to GA, allowing the DELLA protein to bind to the gibberellin responsive gene and inhibit its activity, thereby inhibiting plant growth. When GA is at a high level, GID1 can sense the GA signal, forming GID1-GA-DELLA to degrade DELLA protein, which inhibits the repressing of DELLA on GA signaling [89]. In Q2, the co-expression network of GA signal transduction contained 8 matched lncRNA-mRNA pairs, including 8 lncRNAs and 4 mRNAs. In Qinyou8, the co-expression network of GA signal transduction contained 31 matched lncRNA-mRNA pairs, which included 30 lncRNAs and 4 mRNAs. Two mRNAs (BnaA07g19530D and BnaCnng55170D) co-expressed with lncRNAs, which were down-regulated in both genotypes and which respond to drought stress and re-watering, were annotated to GID1. Downregulated GID1 genes prevented the formation of complexes with GA and DELLA proteins, resulting in the binding of the DELLA protein to the gibberellin response gene, thereby inhibiting seedling growth.
Abscisic acid (ABA) as a signal molecule for plants to perceive stress [90], plays an important role in preventing plant water loss, regulating stomatal opening, and maintaining the balance of cell permeability [90]. ABA binds its receptor PYR/PYL/RCAR (pyrabactin resistant/ PYR-like/regulatory component of ABA) and inhibits the activity of PP2C (protein phosphatases type-2C), which leads to the autophosphorylation of downstream SnRK2 (sucrose non-fermenting 1-related subfamily 2 kinases) and the phosphorylation of downstream ABF transcription factors, regulating the expression of stressrelated genes [91,92]. BnaC07g44670D is homologous to gene ABF (AT4G34000) in Arabidopsis thaliana, which has been reported to be an important gene involved in ABA signaling [93]. In Q2, the co-expression network of ABA signal transduction contained 119 matched lncRNA-mRNA pairs, including 37 lncRNAs and 24 mRNAs. In Qinyou8, the co-expression network of ABA signal transduction contained 207 matched lncRNA-mRNA pairs, which included 73 lncRNAs and 25 mRNAs. In our research, we identified that lncRNAs that co-expressed with BnaC07g44670D, differed between the two genotypes. XLOC_074677, XLOC_ 093758, XLOC_044363 and XLOC_076449, which coexpressed with BnaC07g44670D, were down-regulated in the two genotypes. XLOC_081156 which co-expressed with BnaC07g44670D, was only down-regulated in Qinyou8. These findings suggest that altered lncRNAs may be involved in "plant hormone signal transduction" and regulated differently in the two genotypes. The up-regulation of ABF in response to drought stress can trigger stomatal closure and seed dormancy [94]. The downregulation of ABF in response to re-watering led to weakened ABA signal, which may alleviate rapeseed seedling growth inhibition by ABA.
Studies have shown that MYB was involved in response to abiotic stress, which could be induced by ABA, to participate in the regulation of waxy synthesis pathway of drought stress response [110], and that it promoted the drought resistance of plants by promoting stomatal closure and reducing leaf water loss [111,112]. At present, the research on the possible role of bHLH TFs in plant response to drought stress mainly focuses on stomatal development, trichome development, root hair development, and abscisic acid (ABA) sensitivity [99]. The bHLH-type transcription factor AtAIB depended on ABA signal transduction pathway to participate in the drought resistance response in Arabidopsis [113]. It was found that overexpression of OsbHLH148 in rice induced up-regulation of OsDREB, OsJAZ and other related genes involved in stress response, and in the jasmonic acid signaling pathway, indicating that OsbHLH148 regulated the expression of jasmonic acid signaling pathway-related genes as a starting response factor during drought stress [114]. Among expressed TFs, the most specifically expressed in Q2 and Qinyou8 were MYB and bHLH, respectively. It may be one of the important reasons for the different regulation modes of the two genotypes' response to drought stress and re-watering. Nuclear factor Y (NF-Y) is composed of three distinct subunits (NF-YA, NF-YB, and NF-YC). We found that the Arabidopsis thaliana NFYA5 transcript is strongly induced by abscisic acid (ABA)dependent manner under drought stress, and, the overexpressing of NFYA5 in Arabidopsis thaliana resisted drought stress by controlling stomatal aperture so as to reduce leaf water loss [115]. In this study, NFY accounted for the largest proportion of co-expressed TFs in the two genotypes, respectively. In summary, the two genotypes have different ways of responding to drought stress and re-watering, which is conducive to understanding the molecular regulatory mechanism in response to drought stress, and strengthening our understanding of drought regulatory network.
LncRNA HID1 (HIDDEN TREASURE 1) has been proved to be an important participant in seeding light morphogenesis by regulating PIF3 (phytochrome-interacting factor 3) expression [116]. In Chinese cabbage (Brassica rapa ssp. chinensis), some TFs were cisregulated by the response of lncRNAs to heat stress [47]. Under water stress and during recovery, 189 TFs corresponded to 163 differentially expressed lncRNAs in C. songorica, and there was a bZIP gene predicted to be the target gene of an lncRNA (MSTRG.17203.1) [72]. These studies indicated that there was a regulatory relationship between lncRNAs and TFs. In total, 57 and 94 TFs related to 20 and 24 different families showed coexpression with lncRNAs in two genotypes, respectively. Though, the number of TFs and TF families coexpressed in Qinyou8 higher than Q2, but the occurrence pattern was comparable. The TFs related to the HSF, NF-YA, ERF, bHLH, MYB, GATA, and bZIP families were highly represented in Q2. Similarly, HSF, NF-YA, ERF, bHLH, MYB, WRKY, and bZIP TF families were more enriched in Qinyou8. Among specifically expressed TFs in Q2, a PAT1 gene (BnaC07g49170D) was predicted to be XLOC_096112 target gene and a TGA3 (BnaC05g17700D) was predicted to be XLOC_ 032712 target gene. In Qinyou8, a bHLH69 gene (BnaC01g07430D) was predicted to be a target gene for 10 lncRNAs. In our research, we also found that a bZIP gene (BnaA09g03330D) was predicted to be a target gene for 7 lncRNAs in Q2 and two bZIP genes (BnaA09g03330D and BnaA09g19470D) were predicted to be the target genes for 35 lncRNAs in Qinyou8. This result suggested that the regulation of lncRNAs might play crucial roles in response to drought stress. This would be the next step to explore.
Other lncRNAs involved in drought stress and re-watering Some other candidate functional and regulatory lncRNAs have been detected in response to drought stress and re-watering. We identified that XLOC_052298 and XLOC_094954 were down-regulated in the tolerant genotype and up-regulated in the susceptible genotype, XLOC_012868 was up-regulated in the tolerant genotype and down-regulated in the susceptible genotype. It was noticed that some mRNAs which were co-located with three lncRNAs, were mainly categorized into two categories, i.e. signal transport and defense/stress response.
Drought signals may be perceived by changes in membrane receptor activity. At this time, extracellular signals are converted into intracellular signals, which can lead to the production of second messengers such as Ca 2+ , sugars, ROS and IP 3 delivery systems [117], triggering phosphorylation/dephosphorylation reactions and transmitting information, thereby activating specific transcription factors. After binding to the corresponding cis-acting elements, transcription factors regulate the expression of drought-stress-responsive genes [118]. Serine/threonine protein phosphatase is one of the major enzymes that catalyze the dephosphorylation of proteins [119]. A previous study has demonstrated that serine/threonine protein phosphatase is related to the regulation of anti-reverse signal transduction induced by abscisic acid in plants [120,121]. As the core component of BR signaling, the BES1/ BZR1 transcription factors are activated by the BR signal, bind to the E-box (CANNTG) or BRRE element (CGTG T/CG) of the growth and development-related genes promoter and regulate target gene expression [122][123][124]. BRs, an important plant hormone, improves drought resistance of plants by improving plant osmotic regulation and influencing the activities of antioxidative enzymes [125,126]. Under drought stress, the accumulation of soluble sugars, such as trehalose, has the function of stabilizing the proteins and cell membranes, which is beneficial for the regulation of the balance between the osmotic pressure and the outside of the plant cells [127,128]. Plants with reduced gibberellin (GA) activity, and therefore reduced transpiration, suffer less from leaf desiccation, thereby maintaining higher capabilities and recovery rates [129]. In this study, BnaC02g25020D, BnaC02g25150D and BnaC02g25200D, which co-locate with XLOC_052298, were associated with alpha-trehalose-phosphate synthase, peroxidase, and the BES1/BZR1 homolog protein, respectively. BnaC09g24140D, which co-locates with XLOC_ 094954, was associated with serine/threonine-protein phosphatase. BnaA03g47140D and BnaA03g47400D, which co-locate with XLOC_012868, were associated with superoxide dismutase, gibberellin oxidase, respectively. BnaA03g47370D and BnaA03g47380D, which co-locate with XLOC_012868, were associated with bHLH. Therefore, we believe that these lncRNAs may be related to drought stress and re-watering. However, our knowledge about the potential functions of these dysregulated lncRNAs in response to drought, remains limited. Thus, further investigation is of. great importance.

Conclusion
In this study, 5546 down-regulated and 6997 upregulated mRNAs were detected in Q2, as compared to 7824 and 10,251 in Qinyou8, respectively; 369 downregulated and 108 up-regulated lncRNAs were detected in Q2, compared with 449 and 257 in Qinyou8, respectively. In addition, the interaction networks between lncRNAs and mRNAs were constructed and the function of lncRNAs was then investigated based on the lncRNA-mRNA in-teraction networks. This study found that 4 lncRNAs were annotated significantly to the ABA signaling pathway within a KEGG pathway "plant hormone signal transduction", in Q2, under drought stress and rewatering. Eight mRNAs, which co-locate with three lncRNAs, were mainly categorized into signal transport and defense/stress response under drought stress and re-watering. At the same time, the photosynthesisassociated genes were commonly up-regulated by drought stress and re-watering treatment in Qinyou8. In conclusion, the foregoing outcome indicates that drought stress and re-watering affects the expression of some lncRNAs, and the inter-regulation of lncRNAs and mRNAs may elicit response to drought stress and rewatering. While these findings provide newfound information regarding the potential role of lncRNAs in response to drought stress and re-watering, further research is required to elucidate the molecular mechanisms of significantly dysregulated lncRNAs. The coexpression network suggests that the inter-regulation of lncRNAs and mRNAs is involved in responses to drought stress and re-watering.

Plant materials, growth conditions, and treatments
Seeds of two contrasting rapeseed (Brassica napus L.) genotypes, Q2 (drought-tolerant) and Qinyou8 (droughtsensitive) were obtained from Oil Crops Research Institute (OCRI), Chinese Academy of Agricultural Sciences (CAAS), Wuhan, China. These two contrasting rapeseed genotypes were selected by analyzing the photosynthetic rate, chlorophyll content, carotenoid content, malondialdehyde content, and antioxidants activity of leaves under water stress [130,131]. Under drought stress, Q2 had a relatively higher net photosynthetic rate, the relative water content (RWC), chlorophyll content, carotenoid content, and antioxidants activity in leaves than Qinyou8 [130][131][132][133].
The experiment was conducted in a greenhouse at 25°C, with a photoperiod of 16 h of light and 8 h of darkness, in June 2017, and a humidity rate of 83%. The detailed preparation of the seeds and soils in pots were according to in Xiong et al. [132]. All pots were watered to 75% FC for 18 d (the three-leaf stage) with daily watering, before being subjected to drought stress. Experiment treatment conditions were as follows: (1) 18 days old plants were subjected to water deficit by leaving them un-watered for 8 days (set as drought stress, DS); (2) 18 days old plants were subjected to water deficit by leaving them un-watered for 7 days down to 35% FC [133] and then re-watered for 1 day to 75% FC (set as re-watering, RW). The experiment was carried out using a completely randomized design with three replications. After 8 d of treatment, the 3rd leaves were separately sampled from 5 individuals under each treatment from each replicate (12 samples in total) and quickly stored individually in liquid N 2 .

Determination of physiological parameters
The uniform seedlings of each replicate under the WW, DS and RW treatments were used to measure fresh weight.

RNA extraction, library construction, and Illumina sequencing
The process of RNA extraction and purity were according to Hu et al. [94]. A amount of 3 μg RNA per sample were used to generate cDNA libraries and sequenced. The qualified cDNA libraries were constructed by PCR enrichment and sequenced on a HiSeq X Ten with a sequencing read length of PE125. The 12 gene expression libraries were named DSQ2-1, DSQ2-2, DSQ2-3, RWQ2-1, RWQ2-2, RWQ2-3, DSQinyou8-1, DSQi-nyou8-2, DSQinyou8-3, RWQingyou8-1, RWQinyou8-2, and RWQinyu8-3. The library preparation and deep sequencing were performed by the Novogene Bioinformatics Technology Cooperation (Beijing, China). All the clean reads, obtained after the quality-control process, were deposited in the NCBI Sequence Read Archive with the ID PRJNA574049 for data analysis, as given in the following section.

Mapping to the reference genome
Reference genome and gene model annotation files were downloaded from a genome website (http://brassicadb. org/brad/datasets/pub/Genomes/) directly. Index of the reference genome was built using bowtie v2.0.6 and paired-end clean reads were aligned to the reference genome using TopHat v2.0.9

Quantification of gene expression level
Cuffdiff was used to calculate FPKMs (fragments per kilo-base of exon per million fragments mapped) of both lncRNAs and coding genes in each sample [134]. Values were calculated based on the length of the fragments and reads count mapped to this fragment. Gene FPKMs were computed by summing the FPKMs of the transcripts in each gene group.

Differential expression analysis
Cuffdiff provides statistical routines for determining differential expression in digital transcript or gene expression data using a model based on the negative binomial distribution [134]. The expression strength of each gene was measured by the FPKM method [134] and calculated through averaging expression data of three replicates. The differentially expressed lncRNAs and mRNAs between samples were confirmed by Cufflinks software with the shreshold of q value ≤0.05 and a |log2(FPKM) ratio| ≥ 1. The calculation of q value was according to Trapnell et al. [134]. The treatments of RWQ2/DSQ2 and RWQinyou8/DSQinyou8 were performed.

Construction of the lncRNA-mRNA co-expression network
LncRNA-mRNA co-expression networks were constructed to identify the interactions between proteincoding genes and lncRNAs according to the normalized signal intensities of the specific expression in genes and lncRNAs [135]. We constructed the lncRNA-mRNA coexpression network according to   [136]. Firstly, the expression values of the differentially expressed lncRNAs and mRNAs were obtained. Secondly, the correlation between the differentially expressed lncRNAs and mRNAs was evaluated using the Pearson's correlation coefficient (PCC) from matched mRNA and lncRNA expression profile data. The lncRNA-mRNA pairs with |PCC value ≥0.95| and p < 0.05 were selected as co-regulated lncRNA-mRNA pairs. Subsequently, the network was constructed, in which nodes were lncRNAs or mRNAs. In total, the lncRNA-mRNA co-expression networks were initially constructed based on co-expressed lncRNA-mRNA pairs in each comparison (RWQ2/DSQ2, RWQinyou8/DSQi-nyou8). Ultimately, to visually display the relationship between lncRNAs and target protein-coding RNAs, the interactive networks were constructed using Cytoscape software (3.7.1), (an open source software platform for visualizing complex networks available from http:// cytoscape.org/).

Function classification of the target genes of differentially expressed lncRNAs
Gene Ontology (GO) enrichment analysis of the target genes which co-expressed with differentially expressed lncRNAs were implemented using the GOseq R package, in which gene length bias was corrected. GO terms with corrected p value less than 0.05 were considered significantly enriched with differential expressed genes. We used KOBAS software to test the statistical enrichment of the target genes which co-expressed with differentially expressed lncRNAs in KEGG pathways. The most enriched KEGG was enlisted in order according to the corrected p value. A corrected p value < 0.05 was required for differences to be considered statistically significant.

qRT-PCR analysis
After treated with RNase-free DNase, the RNA samples were used to generate cDNA by using the RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA). Real-time PCR was performed on the ABI 700 platform with the SYBR Green PCR Master Mix system (Takara Co. Ltd., Japan). The 10 μl reaction volume in each well contained 0.5 ng cDNA, 2.5 μl of a mixture containing 1.2 μM each of the forward and reverse primers and 5 μl of master mix. The PCR amplification procedures were as: one iniative cycle of 30 s at 95°C; followed by denaturation at 94°C for 30 s, primer annealing at 60°C for 30 s, and then extension at 72°C for 1 min; finally, an extra extension at 72°C for 10 min. The primer sequences for the randomly selected lncRNAs were shown in Additional file 9. Each PCR reaction were repeated three times independent and the expression strength of each lncRNA was set as their average value.