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 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 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-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 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 [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 re-watering 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 co-expression 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 BnaA05g14780D (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, up-regulated XLOC_017878, targeting up-regulated BnaA09g42140D 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 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 [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 BnaC06g18770D 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 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 [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. 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 [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 stress-related 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 co-expressed 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 down-regulation of ABF in response to re-watering led to weakened ABA signal, which may alleviate rapeseed seedling growth inhibition by ABA.
Transcription factors functioned under drought stress and re-watering
Transcription factors have been confirmed to play a crucial role in regulating drought stress in plants [89, 95, 96]. Previously, MYB [97, 98], bHLH [99 100], WRKY [101], ERF [102], NFY [103], GATA [104], PIF [105, 106], ABA-INDUCIBLE BHLH-TYPE TRANSCRIPTION FACTOR (AIB) [107], HSF [108], and bZIP [109] had been proposed as being responsive to abiotic stresses. In this study, these TFs were induced to express under drought stress and re-watering.
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 cis-regulated 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 co-expression with lncRNAs in two genotypes, respectively. Though, the number of TFs and TF families co-expressed 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 Ca2+, sugars, ROS and IP3 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 (CGTGT/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.