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


 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) in response to drought stress and rehydration treatment at the seedling stage.

Results: A total of 5,546 down-regulated and 6,997 up-regulated mRNAs were detected in Q2 compared with 7,824 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. We further identified 34 TFs corresponding to 126 differently expressed lncRNAs in Q2, and 45 TFs corresponding to 359 differently 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 and co-located with BnaC07g44670D, which are associated with plant hormone signal transduction. Additionally, some mRNAs which were co-located with LNC_002535 (XLOC_052298), LNC_004924 (XLOC_094954) and LNC_000539 (XLOC_012868) were mainly categorized as signal transport and defense/stress response. Finally, co-expression network analysis indicated that the co-expression network of Q2 was composed of 145 network nodes and 5,175 connections, while the co-expression network of Qinyou8 was composed of 305 network nodes and 22,327 connections.

Conclusions: The differentially expressed mRNAs and lncRNAs may play important roles in response to drought stress and rehydration treatments and could provide basic information for new ways to improve the drought resistance of Rapeseed Brassica napus .

5 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 are still lacking, especially in B. napus. In order to further understand the molecular mechanisms of the response of B. napus to drought-stress, 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 determined whether mRNAs and lncRNAs are expressed differentially in plants subject to drought stress and rehydration treatments. Therefore, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to predict the biological roles and potential signaling pathways of these differentially expressed lncRNAs. Moreover, lncRNA-mRNA network analysis was conducted to further explore the potential roles of differentially expressed lncRNAs in response to drought stress.

Differently expressed lncRNAs and mRNAs under drought stress and re-watering
In this study, we extracted RNAs from 12 samples (two treatments, two test materials, three biological replicates) and tested their quality before performing RNA sequencing (Table S1). 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 (Table S2). 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 qvalue < 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 6 down-regulated, 257 up-regulated) of Qinyou8 were specifically dysregulated after rewatering ( Fig. 1). In addition, there were 12,543 mRNAs (5,546 down-regulated, 6,997 upregulated) and 18,075 mRNAs (7,824 down-regulated, 10,251 up-regulated) differentially expressed in Q2 and Qinyou8 (Fig.1), respectively.

qRT-PCR Validation
To validate our results independently and determine the role of lncRNAs in response to drought stress, we selected nine lncRNAs which were differentially expressed in two genotypes. As shown in Fig. 2, the linear relationship between the expression profiles of lncRNAs, is shown by qRT-PCR and RNA-seq, and the line is fitted using linear regression (R 2 = 0.91519, slope = 0.91646). The real-time PCR results verify the expression patterns obtained with transcriptome sequencing, indicating that the lncRNAs expression based on RNA-seq data are reliable.

GO analysis of the differently expressed lncRNAs and mRNA under drought stress and re-watering
Studies have shown that lncRNA can indirectly affect the expression of mRNA, and can also directly bind to mRNA, thus affecting translation [47,48,49], shearing [50,51], and degradation of mRNA [52]. 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 the mRNAs co-expressed with the differentially expressed lncRNAs. This analysis was performed to determine the major molecular functions, biological processes, and cellular components with which the differentially expressed lncRNAs were associated.
The down-regulated mRNAs, that were co-expressed with lncRNAs, were assigned to 32 and 34 significant terms in Q2 and Qinyou8 (Fig. 3a), respectively. For down-regulated mRNAs co-expressed with lncRNAs in Q2, the most significant GO terms for biological 7 process were oxidation-reduction process (GO:0055114), protein dephosphorylation (GO:0006470), dephosphorylation (GO:0016311), 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 lncRNAs, the important GO terms for biological process were single-organism metabolic process (GO:0044710), oxidationreduction 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 CCAATbinding factor complex (GO:0016602) were the most significantly enriched GO terms.
The up-regulated mRNAs co-expressed with lncRNAs were assigned to 23 and 31 significant GO terms in Q2 and Qinyou8 (Fig. 3b), respectively. Of the enriched GO terms in the biological process category for up-regulated mRNAs co-expressed with lncRNAs in Q2, primary metabolic process (GO:0044238) was the most dominant group, followed by 8 organic substance metabolic process (GO:0071704), macromolecule metabolic process 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 stresses.

Identification of transcription factors under drought stress and re-watering
Under drought stress, transcription factors (TFs) can be used as regulators of drought stress, and they would bind to cis-acting elements in the promoter region of related genes to regulate the expression of downstream genes [53]. In our research we found that there were 334 differentially expressed genes in Q2, and 487 differentially expressed genes in Qinyou8; when compared to the transcription factor database of Arabidopsis, 211 TFs were found to be co-expressed in two genotypes, and 12 TFs were conversely expressed in two genotypes. In addition, 123 TFs that were specifically expressed in Q2 were classified into 38 groups; 10 of these TF families comprised 69.92% of these groups, including MYB (23 TFs), basic helix-loop-helix (bHLH) (18 TFs), ERF (11 TFs), WRKY (8 TFs), PIF (7 TFs), GATA Additionally, 6 groups conversely expressed between two genotypes (con-expression TFs); these included bHLH (4 TFs), NFY (4 TFs), DIVARICATA (1 TFs), ERF (1 TFs), PIF (1 TFs) and TCP (1 TFs) (Fig. 5d). We also analyzed the relationship between mRNAs co-expressed with lncRNAs and transcription factors. In Q2, 34 TFs corresponded to 126 differently expressed lncRNAs (Table S3), while there were 45 TFs corresponding to 359 differently expressed lncRNAs in Qinyou8 (Table S4).

LncRNA-mRNA Co-Expression Network
To further characterize the role of differentially expressed lncRNAs, we identified the network of lncRNAs by co-expression analysis. Co-expression network analysis indicated that the co-expression network of Q2 was composed of 145 network nodes and 5,175 connections, while the co-expression network of Qinyou8 was composed of 305 network nodes and 22,327 connections.
One lncRNA may regulate multiple protein-coding genes, and vice versa [54,55,56]. 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 coexpression network of Qinyou8 indicated that 1 mRNA may correlate with 1 to 43 lncRNAs, and 1 lncRNA may correlate with 1 to 375 mRNAs. LNC_003598 (XLOC_071559) was the largest node in the network in both genotypes, respectively.

Discussion
Several recent studies have revealed that lncRNAs play an important role in response to drought stresses [31,39,41,57]. However, a genome-wide identification and characterization of lncRNAs responses to drought stress and re-watering treatments are still lacking, especially in B. napus. Accordingly, we constructed lncRNA and mRNA libraries, and annotated, identified, and verified those lncRNAs that are involved in drought stresses.

Differential mRNAs and lncRNAs expression in two contrasting genotypes under drought stress and re-watering
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, 5,546 down-regulated and 6,997 up-regulated mRNAs were detected in Q2 compared to 7,824 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 up-regulated) in both genotypes, among which, 1 lncRNA (LNC_000539, 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. 6). 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. Additionally, 9 lncRNAs that we had identified 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,58], response to stimulus [59], response to abiotic stimulus [58], biosynthetic process [58], structural constituent of ribosome [56], photosynthesis [60] and oxidoreductase activity [60], 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 [61], and plant hormone signal transduction [42,62]. 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-13 expressed with 40 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 up-regulated mRNAs co-expressed with lncRNAs belonging to photosynthesis were significantly and specially enriched in Qinyou8.
A total of 7 and 5 mRNAs co-expressed with lncRNAs were assigned into photosynthesis and photosynthesis-antenna proteins, respectively. The genes involved in photosynthesis were generally down-regulated by drought [63,64], our results indicate that the photosynthesis-associated genes were commonly up-regulated by drought stress and rewatering treatment. These results indicate that lncRNAs could play a role in many biological processes responding to drought stress 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.
To further explore the different regulation mechanisms between the two genotypes in response to drought stress and re-watering, we focused on the analysis of plant hormone signal transduction. In this pathway, there were 36 mRNAs co-expressed with 40 lncRNAs in Q2 and 51 mRNAs co-expressed with 119 lncRNAs in Qinyou8, respectively. The differentially expressed mRNAs, that were co-expressed with lncRNAs, were involved in auxin, cytokinin, gibberellin, and abscisic acid signaling pathways in both the genotypes.
Some mRNAs related to the ethylene and salicylic acid signaling pathways were specifically expressed in Q2, while mRNAs involved in the two signaling pathways of brassinosteroid and jasmonic acid were specifically expressed in Qinyou8. Among these signaling pathways, more of the mRNAs, which co-expressed with lncRNAs, were associated with the ABA signaling pathways than those of other phytohormones, which is 14 in good agreement with previous studies that had considered ABA to be an early warning signal for plant drought responses [65,66].
Auxin (IAA) as a phytohormone, is essential for signaling, transport, growth and development of a plant [67]. 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 [68]. 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 [69]. 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 [70]. The GH3 gene encodes an auxin-binding enzyme that acts as a feedback regulator of auxin by reducing the level of beneficial auxin [71]. Studies have shown that SAUR has a positive regulatory effect on cell growth [72]. 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), GH3 (1 differentially expressed mRNA co-expressed with lncRNAs in Q2, and 3 differentially expressed mRNAs co-expressed with lncRNAs in Qinyou8) and SAUR (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, which may promote vegetative growth through cell enlargement.

Cytokinin (CK) plays an important role in various physiological functions in plants, such
as promoting cell division, inducing shoot formation and promoting its growth [67].
Cytokinin signaling is based on a two-component 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 [73]. 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 [74]. One gene (BnaA01g17750D co-expressed with lncRNAs) was annotated to B-ARR and down-regulated in both the genotypes, while only 1 gene (BnaC06g18770D co-expressed with lncRNAs) was annotated to A-ARR in Qinyou8. These indicate that rehydration after drought is beneficial to the growth of seedlings.
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 [75]. GIBBERELLIN INSENSITIVE 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 [76]. 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 [77]. 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 [76]. 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. Down-regulated 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 [78], plays an important role in preventing plant water loss, regulating stomatal opening, and maintaining the balance of cell permeability [79]. 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 stress-related genes [80,81]. BnaC07g44670D is homologous to gene ABF (AT4G34000) in Arabidopsis thaliana, which has been reported to be an important gene involved in ABA signaling [82]. In our research, we identified that lncRNAs that co-expressed and were colocalized with BnaC07g44670D, differed between the two genotypes. LNC_003769 (XLOC_074677), LNC_004839 (XLOC_093758), LNC_002030 (XLOC_044363) and LNC_003865 (XLOC_076449), which co-expressed with BnaC07g44670D, were downregulated in the two genotypes. LNC_004119 (XLOC_081156), which co-expressed with BnaC07g44670D, was only down-regulated in Qinyou8, while LNC_004200 (XLOC_082728) and LNC_004199 (XLOC_082727), which co-localized with BnaC07g44670D, were only upregulated in Qinyou8. These findings suggest that altered lncRNAs may be involved in "plant hormone signal transduction" and regulated differently in the two genotypes. The down-regulation of ABF in response to drought stress and re-watering can trigger stomatal opening and promote normal physiology.
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 [96], and that it promoted the drought resistance of plants by promoting stomatal closure and reducing leaf water loss [97,98]. 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 [87]. The bHLH-type transcription factor AtAIB depended on ABA signal transduction pathway to participate in the drought resistance response in Arabidopsis [99]. 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 [100]. 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 [101]. 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, 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 [102]. In Chinese cabbage (Brassica rapa ssp. chinensis), some TFs were cisregulated by the response of lncRNAs to heat stress [103]. 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) [60]. These studies indicated that there was a regulatory relationship between lncRNAs and TFs. We found that there were 34 TFs corresponding to 126 differently expressed lncRNAs in Q2, while there were 45 TFs corresponding to 359 differently expressed lncRNAs in Qinyou8. Among specifically expressed TFs in Q2, a PAT1 gene (BnaC07g49170D) was predicted to be LNC_004990 (XLOC_096112) target gene and a TGA3 (BnaC05g17700D) was predicted to be LNC_001409 (XLOC_032712) target gene. In Qinyou8, a bHLH69 gene (BnaC01g07430D) was predicted to be a target gene for 10 lncRNAs. 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 LNC_002535 (XLOC_052298) and LNC_004924 (XLOC_094954) were down-regulated in the tolerant genotype and up-regulated in the susceptible genotype, LNC_000539 (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 [104], 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 [105]. Serine/threonine protein phosphatase is one of the major enzymes that catalyzes the dephosphorylation of proteins [106]. 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 [107,108]. 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 (CGTGT/CG) of the growth and development-related genes promoter and regulate target gene expression [109,110,111]. Brassinosteroid (BRs), as an important plant hormone, improves drought resistance of plants by improving plant osmotic regulation and influencing the activities of antioxidative enzymes [112,113]. 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 [114,115]. Plants with reduced gibberellin (GA) activity, and therefore reduced transpiration, suffer less from leaf desiccation, thereby maintaining 20 higher capabilities and recovery rates [116]. In this study, BnaC02g25020D, BnaC02g25150D and BnaC02g25200D, which co-locate with LNC_002535, were associated with alpha-trehalose-phosphate synthase, peroxidase, and the BES1/BZR1 homolog protein, respectively. BnaC09g24140D, which co-locates with LNC_004924, was associated with serine/threonine-protein phosphatase. BnaA03g47140D and BnaA03g47400D, which co-locate with LNC_000539, were associated with superoxide dismutase, gibberellin oxidase, respectively. BnaA03g47370D and BnaA03g47380D, which co-locate with LNC_000539, 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.

LncRNAs and mRNAs are involved in an interaction network
In our study, an lncRNA-mRNA network analysis was employed to identify interactions between differentially expressed lncRNAs and differentially expressed mRNAs, as previously described [55,88]. In Q2, 40 lncRNAs interacted with 36 mRNAs in plant hormone signal transduction (Fig.7a); 15 lncRNAs interacted with 10 mRNAs in the GO of the response to abiotic stimulus (Fig.7b); and 8 lncRNAs interacted with 6 mRNAs in the GO of response to water stimulus (Fig. 7c). In Qinyou8, 5 lncRNAs interacted with 7 mRNAs in photosynthesis (Fig. 7d); 68 lncRNAs interacted with 16 mRNAs in the GO of response to abiotic stimulus (Fig. 7e); and, 58 lncRNAs interacted with 11 mRNAs in the GO of response to water stimulus (Fig. 7f). The response to drought stress in plants is a complicated process, involving several genes and metabolic networks; hormones synthesis is one of the important factors involved [39]. Since these processes (response to abiotic stimulus, response to water stimulus) are important to regulate plant adaptation to drought stress, and the lncRNAs potentially regulate these genes, it could be proposed 21 that lncRNAs play important roles in response to drought stress. The co-expression network suggests that the inter-regulation of lncRNAs and mRNAs is involved in the response to drought and warrants further study.

Conclusions
In this study, 5,546 down-regulated and 6,997 up-regulated mRNAs were detected in Q2, as compared to 7,824 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. 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 re-watering. 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 photosynthesis-associated 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 re-watering. While these findings provide newfound information regarding the potential role of lncRNAs in response to drought stress and rewatering, further research is required to elucidate the molecular mechanisms of significantly dysregulated lncRNAs. The co-expression network suggests that the interregulation of lncRNAs and mRNAs is involved in responses to drought stress and rewatering.

Plant materials and experimental conditions
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%, at the Experimental Station of the Oil Crops Research Institute (OCRI), Chinese Academy of Agricultural Sciences (CAAS), Wuhan, China. The seeds of the two B. napus genotypes, Q2 (droughttolerant) and Qinyou8 (drought-sensitive) were obtained from Oil Crops Research Institute (OCRI), Chinese Academy of Agricultural Sciences (CAAS), Wuhan, China. Their resistance to drought were screened and confirmed [95,[117][118][119][120]. The plant material used in our study has not been deposited in a publicly available herbarium. The seeds were surfacesterilized with 0.1% HgCl 2 for 600 s and rinsed with sterile running water for 1200 s. Ten seeds were then sown in plastic pots (height 25 cm, width 15 cm) that were filled with 1.0 kg of soil, vermiculite and sand mixture (soil: vermiculite : sand = 2:1:1, v/v/v). The vermiculite mixture had a field capacity (FC) of 45.21%, FC being measured according to Wilcox (1965) [121] and Duan et al. [122]. Four plants were placed at a uniform distance 10 days after sowing. Each pot was supplemented with 5 g of fertilizer containing N-13%, P-10%, and K-14% before sowing, and was weighed daily to calculate the soil water content.
All pots were watered to 75% FC for 18 d (the three-leaf stage) with daily watering, before being subjected to drought stress. Then, two water regimes were imposed: in the DS treatment, water was withheld, while in the RW treatment, the plants were re-watered to 75% FC after drought stress. Two treatments were conducted as follows: (1) DS-natural drought to 35% FC (about 3 days) [118], continued drought for 8 days; (2) RW-natural drought to 35% FC (about 3 days), continued drought for 7 days, and then re-watered for 1 day to 75% FC. The experiment was carried out using a completely randomized design with three replications.

RNA preparation
Total RNA was extracted from the 3 rd leaves of each sample using TRIzol reagent (Invitrogen, Burlington, ON, Canada) according to the manufacturer's instructions. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the Nanophotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). In the reaction buffer, dNTPs with dTTP were replaced by dUTP. Remaining overhangs were converted into blunt ends through exonuclease/polymerase activities. After adenylation of 3' ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150-200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB USA) was used on size-selected, adaptorligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index 24 (X) Primer. Finally, the products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) [123].

Library preparation for lncRNA sequencing
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 [124]. 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 [124]. Transcripts with a P-adjust < 0.05 were assigned as differentially expressed. Differentially expressed lncRNAs and mRNAs were identified through the random variance model with P values calculated by the paired t test. The significance cutoff for the up-regulated and down-regulated lncRNAs and mRNAs was a fold change ≥ 2.0 with P ≤ 0.05.

GO and KEGG enrichment analysis
Gene Ontology (GO) enrichment analysis of differentially expressed genes or lncRNA target genes 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 differential expression genes or lncRNA target genes in KEGG pathways.

LncRNA-mRNA co-expression network
Pearson correlation coefficient (PCC) was used to depict the co-expression relationship between lncRNA and mRNA according to their expression levels. lncRNA-mRNA pairs with |PCC value ≥ 0.90| and p < 0.05 were retained for network construction, which was deciphered using Cytoscape 3.7.1, (an open source software platform for visualizing complex networks available from http://cytoscape.org/).

qRT-PCR analysis
RNA samples were treated with RNase-free DNase and reverse-transcribed using the RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA). Real-time PCR was carried out using the SYBR Green PCR Master Mix system (Takara Co. Ltd., Japan) on an ABI 700 Real-time PCR platform. Reactions (10 μl) were carried out in 96-well format and contained approximately 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 conditions were set as follow: one cycle of 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, and 60°C for 30 s. Primers for PCR assay that were used for candidate genes were designed using the Primer Premier 5 software. The primer sequences for the selected lncRNAs are listed in Table S5. All quantitative real-time PCR reactions for each lncRNA were performed using three independent biological replicates, and each was based on three technical repeats.

Consent for publication
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Table S1. RNA-seq data for 12 samples. Table S2. RNA amount obtained from each treatment for RNA-seq analysis. Table S3. Transcription factors corresponding to differently expressed lncRNAs in Q2. Table S4. Transcription factors corresponding to differently expressed lncRNAs in Qinyou8. Table S5. List of primers for quantitative real-time PCR.

Supplementary Material
44 Figure 1 The numbers of differentially expressed lncRNAs and mRNAs in two genotypes (Q2 and Qinyou 8) in response to drought stress and re-watering treatments.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to