Genome-Wide Differential DNA Methylation and miRNA Expression Profiling Reveals Epigenetic Regulatory Mechanisms Underlying Nitrogen-Limitation-Triggered Adaptation and Use Efficiency Enhancement in Allotetraploid Rapeseed

Improving crop nitrogen (N) limitation adaptation (NLA) is a core approach to enhance N use efficiency (NUE) and reduce N fertilizer application. Rapeseed has a high demand for N nutrients for optimal plant growth and seed production, but it exhibits low NUE. Epigenetic modification, such as DNA methylation and modification from small RNAs, is key to plant adaptive responses to various stresses. However, epigenetic regulatory mechanisms underlying NLA and NUE remain elusive in allotetraploid B. napus. In this study, we identified overaccumulated carbohydrate, and improved primary and lateral roots in rapeseed plants under N limitation, which resulted in decreased plant nitrate concentrations, enhanced root-to-shoot N translocation, and increased NUE. Transcriptomics and RT-qPCR assays revealed that N limitation induced the expression of NRT1.1, NRT1.5, NRT1.7, NRT2.1/NAR2.1, and Gln1;1, and repressed the transcriptional levels of CLCa, NRT1.8, and NIA1. High-resolution whole genome bisulfite sequencing characterized 5094 differentially methylated genes involving ubiquitin-mediated proteolysis, N recycling, and phytohormone metabolism under N limitation. Hypermethylation/hypomethylation in promoter regions or gene bodies of some key N-metabolism genes might be involved in their transcriptional regulation by N limitation. Genome-wide miRNA sequencing identified 224 N limitation-responsive differentially expressed miRNAs regulating leaf development, amino acid metabolism, and plant hormone signal transduction. Furthermore, degradome sequencing and RT-qPCR assays revealed the miR827-NLA pathway regulating limited N-induced leaf senescence as well as the miR171-SCL6 and miR160-ARF17 pathways regulating root growth under N deficiency. Our study provides a comprehensive insight into the epigenetic regulatory mechanisms underlying rapeseed NLA, and it will be helpful for genetic engineering of NUE in crop species through epigenetic modification of some N metabolism-associated genes.


Introduction
Nitrogen (N) is a key macronutrient for plant growth, development, and yield capacity [1]. Nitrate (NO 3 − ) is an N nutrient source that is preferentially absorbed by terrestrial plants [2]. Low nitrate availability is often observed in most agricultural soils; therefore, a huge amount of energy is consumed to annually apply large numbers of chemical N fertilizers to ensure crop production [3]. However, owing to leaching loss, soil surface runoff, denitrification, volatilization, and microbial consumption, more than half of the N nutrients are wasted [4]. Thus, increasing plant adaptability to N limitation is key to reducing N fertilizer application in modern agriculture production, which contributes to development of environmentally friendly agriculture and a sustainable ecosystem. Epigenetic DNA modification, not changing genomic sequences, plays a pivotal role in maintaining genome stability, orchestrating gene expression across plant development, and regulating plant responses to environmental stimulations [5]. Recent studies have highlighted three epigenetic hallmarks: DNA methylation, histone modifications, and modification from small RNAs [6]. In plants, 5-methylcytosine (5 mC) is a major epigenetic signature that is modified by DNA methyltransferases. DNA methylation occurs not only predominantly at symmetric CG sites but also at CHG and asymmetric CHH regions (where H represents A, T, or C) [7]. Nutrient stresses in plants, such as N limitation [8], phosphate starvation [9,10], zinc deficiency [11], and sulfur deficiency [12], have been shown to change global DNA methylation at the whole genome level.
Endogenous microRNAs (miRNAs) are a class of evolutionarily conserved small noncoding RNA molecules that regulate plant transcriptional and post-transcriptional responses to various biotic and abiotic stresses [13]. Some miRNAs modulate plant adaptation to N limitation through modulation of their target genes [14]. Arabidopsis miR169 was downregulated under N deficiency, which induced an accumulation of its targets encoding nuclear factor Y subunit A family members. Plants overexpressing miR169 accumulated less N and showed greater sensitivity to N limitation, which was associated with the downregulation of important nitrate transporter (NRT) genes AtNRT2. 1 and AtNRT1.1 [15]. In addition, the miR167/ARF8 and miR393/AFB3 regulatory circuits were shown to function in N-responsive regulatory networks, which regulated lateral root outgrowth [16]. Repression of an E3 ubiquitin ligase, NLA (Nitrogen Limitation Adaptation), by miRNA827 has been verified at the post-transcriptional (cleavage of mRNA) and translational levels under N deficiency [17]. N starvation triggered a striking decrease in miR164 levels with an increase in the mRNA level of its target, ORE1 [18].
Allotetraploid Brassica napus (A n A n C n C n , 2n = 4x = 38) is an important oleaginous crop worldwide [19]. B. napus is more susceptible to N deficiency than grain crops to maintain optimal growth and development. However, the N use efficiency (NUE) of B. napus is much lower than those of other major crops [20]. Numerous duplicated segments and homoeologous regions are found within the B. napus genome [19], indicating that there might be complicated molecular regulatory networks involving NLA and NUE. To understand and explore epigenetic regulatory mechanisms underlying the adaptive responses of rapeseed plants to N limitation, the aim of our present study is to (i) characterize physiological responses, including carbon-N metabolism, root system architecture (RSA), and ionomic profiling, of rapeseed plants under N limitation, (ii) analyze differentially expressed miRNAs between N sufficiency and N limitation, as well as their targets genes, and (iii) identify differentially methylated regions/genes between N sufficiency and N limitation. Our study will provide helpful information for genetic engineering of NUE and NLA in plants through epigenetic modification of some N-metabolism-associated genes.

Physiological Adaptive Responses of Rapeseed Plants to N Limitation
Under N limitation, rapeseed plants showed smaller and chlorotic leaves ( Figure 1A-D), which was also confirmed by lower biomasses ( Figure 1B). Compared with the decrease in root (R) biomasses, a higher degree of reduction in the shoot (S) biomasses was observed under low N condition ( Figure 1B). Based on the above result, we observed that N limitation led to an increase in the R/S ratio ( Figure 1C). Moreover, limited N, compared with N sufficiency, resulted in over 75% reduction in the total leaf areas ( Figure 1D). Flow of organic acids into amino acids was significantly decreased in plant shoots under N limitation [21]. Under low N, more organic acids were diverted to sugar compounds in the leaves ( Figure 1E-G). Anthocyanin biosynthesis is usually considered as an adaptive response to N limitation [21], and we found that the anthocyanin concentration under N limitation was up to almost twofold of that under N sufficiency ( Figure 1E,F). Moreover, the concentrations of some other carbohydrates, such as glucose, fructose, and sucrose, were also significantly increased under low N ( Figure 1G). Besides great alteration in the shoot performance and metabolism ( Figure 1A-G), N limitation also induced a dramatic change in rapeseed RSA. Although the root biomasses were markedly reduced ( Figure 1B), total root length, root surface area, and root volume were significantly increased under N limitation ( Figure 1H,I). However, root diameters were decreased by more than 40% under low N ( Figure 1I), which might be a key factor causing reduction in root biomasses.
Subsequently, we examined the effect of N limitation on ionomic profiling in rapeseed plants ( Figure 2A). Indeed, low N reduced nitrate concentrations in both shoots and roots ( Figure 2A). However, compared with N sufficiency, a larger proportion of nitrate was translocated from roots to shoots ( Figure 2B), which finally resulted in a significantly higher NUE under N limitation ( Figure 2C). Besides the reduction in nitrate concentrations (Figure 2A), low N induced remarkable decreases in the concentrations of most cations, such as K + , Ca 2+ , Mg 2+ , and Cu 2+ , in both shoots and roots ( Figure 2A). Among them, with exception, the concentrations of Fe 2+ and Mn 2+ in the shoots and Zn 2+ in the roots were not significantly changed between high N and low N conditions (Figure 2A).

Transcriptional Adaptive Responses of Rapeseed Plants to N Limitation
In our previous study [22], we performed a genome-wide transcriptional analysis of rapeseed plants under sufficient N and limited N conditions. We identified a total of 4346 differentially expressed genes (DEGs) in the shoots and 9693 in the roots between the two N supply conditions, respectively [22]. In this study, to better understand molecular responses of rapeseed plants to N limitation, we used the online Panther program to retrieve the functional ontology of identified DEGs. The result showed that transporters (mainly nutrient transporters) and transferases (including glycosyltransferases and methyltransferases) were highly accumulated in both shoots and roots ( Figure 2D). This indicated that nutrient transport and epigenetic modification, such as DNA methylation, might play crucial roles in the adaptive responses of rapeseed plants to N limitation stress.
To examine molecular responses of N-metabolism-associated genes to N deficiency, we determined their differential expression between high N and low N conditions by RT-qPCR assays ( Figure 2E). In terms of N transporter genes that are responsible for N absorption, we focused on the dual-affinity nitrate transporter (NRT) gene, NRT1.1/NPF6.3 [23], and the two-component high-affinity nitrate transport system genes, NRT2.1 and NAR2.1/NRT3.1 [24]. Obviously, the expression of three-type NRTs above-mentioned was significantly induced in the roots after exposure to 72-h N starvation. Furthermore, we investigated the expression of NRTs that are involved in long-distance translocation of N nutrients, such as NRT1.5/NPF7.3 and NRT1.8/NPF7.2, which are responsible for root nitrate xylem loading and unloading [25,26], respectively. The results showed that N depletion elevated the mRNA level of NRT1.5, whereas repressed the NRT1.8 expression. The expression of NRT1.7/NPF2.13, which facilitates N nutrient remobilization from source leaves to sinks [27], was also enhanced under N limitation. Moreover, the expression of a vacuolar nitrate influx transporter gene, CLCa, was lower under N limitation than under N sufficiency. Furthermore, N limitation repressed the expression of the nitrate reductase (NIA/NR) gene NIA1 and increased the mRNA level of the glutamine synthetase (Gln/GS) gene Gln1;1/GS1;1.
Considering that DNA methyltransferase terms were highly enriched in the GO analysis, we decided to investigate differential expression profiling of DNA methyltransferase genes. We identified a total of 12 DNA methyltransferase genes that were differentially expressed in the shoots and roots under N limitation ( Figure 2F). Among the eight DEGs in the shoots, most of them were upregulated, and a similar expression pattern was also observed in the roots ( Figure 2F). Subsequently, we selected a DNA methyltransferase gene (BnaC06g23520D) and examined its expression using RT-qPCR assays; the result confirmed that N limitation led to a significant increase in the DNA methyltransferase expression ( Figure 2G). The finding indicated that N limitation might change DNA methylation status in rapeseed plants.

Genome-Wide High-Resolution DNA Methylation Fingerprints in Response to N Limitation
To further identify the effect of N limitation on genome-wide DNA methylation, we conducted a global comparison of DNA methylation patterns between N sufficiency and N limitation by whole genome bisulfite sequencing (WGBS). The WGBS generated a total of 120 million raw reads for sample library by paired-end sequencing, and more than 90% of them were clean reads, which covered more than 30× depth of the rapeseed genome (Supplementary Table S1). More than 75% cytosines were covered by at least five sequencing reads in the rapeseed genome. Low base calling error rates (Q 20 < 5%, Q 30 < 10%) and high bisulfite conversion rates (> 99%) suggested high quality of the global WGBS data (Supplementary Table S1).
In general, the genome-wide DNA methylation level under the CG context was the highest, followed by that under CHG and CHH contexts ( Figure 3A, Supplementary Figures S1A, S2, and S3). The tendency of variations in DNA methylation across genomic regions was consistent in the shoots and roots under both N sufficiency and N limitation (Supplementary Figure S2). The average methylation levels were uneven among the 19 chromosomes, and the methylation abundance of the A n subgenome was significantly lower than that of the C n subgenome from all the DNA libraries under all the CG, CHH, and CHG contexts ( Figure 3B, Supplementary Figures S1B and S3). To better understand the relationship between DNA methylation and gene expression, we divided the rapeseed genome into several functional categories, including gene promoter, gene body, and 2.0-kb downstream regions of genes. The methylation degrees of gene promoters and downstream regions were obviously higher than those of gene bodies (Supplementary Figure S4), with a slightly higher methylation level for introns than for exons (Supplementary Figure S5). Overview of the genome-wide differential DNA methylation fingerprints in the rapeseed plant shoots between N sufficiency and N limitation. (A) Circos plots showing the density of differentially methylated regions across the rapeseed genome. In the Circos figure, the terms were as follows outside-to-inside: (i) chromosome, (ii) 5mC density in the CG, (iii) CHG, and (iv) CHH contexts, and (v) gene density of each chromosome in alloteraploid rapeseed. (B) Density and level plot of 5mC in CG, CHG, and CHH contexts in the gene bodies on each chromosome. In the Circos figure, the terms were as follows outside-to-inside: (i) chromosome, (ii) 5mC levels in the CG, CHG, and CHH contexts, (iii) gene density, and (iv) 5 mC densities in the CG, CHG, and CHH contexts in the gene bodies. (C-E) Density of hypermethylated/hypomethylated regions in the (C) CG, (D) CHG, and (E) CHH contexts across the rapeseed genome. In the Circos figure, the terms were as follows outside-to-inside: (i) chromosome, (ii) densities of hypermethylated regions in the CG, CHG, and CHH contexts, (iii) gene density, and (iv) densities of hypomethylated regions in the CG, CHG, and CHH contexts. (F-H) Differential methylation levels in the (F) CG, (G) CHG, and (H) CHH contexts across the rapeseed genome. In the Circos figure, the terms were as follows outside-to-inside: (i) chromosome, (ii) hypermethylation levels in the CG, CHG, and CHH contexts, (iii) DNA methylation differences, and (iv) hypermethylation levels in the CG, CHG, and CHH contexts. The times of 0 h and 72 h refer to the times after treatment of N deficiency. To further identify differential DNA methylation profiling between N sufficiency and N limitation, we characterized genome-wide differentially methylated regions (DMRs) in the shoots and roots. From the global view, there were 4474 and 4485 differentially methylated cytosines for CG in the shoots and roots, 3205 and 1200 for CHG, and 18,987 and 4728 for CHH, respectively ( Figure 3C-E, Supplementary Figure S1C-E). In detail, there were high peaks of DMR density on the A2/A9 chromosomes of the A n subgenome and on the C2/C4/C7 chromosomes of the C n subgenome ( Figure 3C-E, Supplementary Figure S1C-E). The majority of DMR lengths under the CG and CHG contexts were fewer than 100 nucleotides (nts), which were obviously shorter than those under the CHH context (Supplementary Figure S6). Different genomic features of differential methylation were observed between N sufficiency and N limitation ( Figure 3F-H and Figure 4, Supplementary Figure  S1F-H). In the CG context, we did not identify significant methylation differences within all the gene body and their 2.0-kb upstream (namely promoter) and downstream regions. In the CHG context, numerous hypermethylation was detected under N limitation compared with that under N sufficiency. However, in the CHH context, a discrepancy in response to N limitation was observed between the shoots and roots; in the shoots, N limitation led to DNA hypomethylation, but the pattern was opposite in the roots. In total, we identified 2389 differentially methylated genes (DMGs), including 1206 hypermethylated genes and 1284 hypomethylated genes, in the shoots between N sufficiency and N limitation. In the roots, we detected 2705 DMGs, including 1511 hypermethylated genes and 1316 hypomethylated genes between N sufficiency and N limitation (Supplementary Figure S7A). The majority of DMGs were identified in the CG genomic context, the number of which was subsequently followed by DMGs in the CHH and CHG regions (Supplementary Figure S7B).
The GO of DMGs and their involvement in KEGG pathways were used to identify the function of DMGs between N sufficiency and N limitation. Cellular and metabolic processes belonging to the biological process item, cell parts and protein complexes belonging to the cellular component group, and transcription binding and transporter activities belonging to the molecular function category were the most overrepresented in the GO analysis ( Figure 5A). Furthermore, we mapped the DMGs onto KEGG pathways in which they were involved ( Figure 5B). In the shoots, the items of starch, sucrose, and anthocyanin (phenylpropanoid pathway) metabolism, photosynthesis, and protein processing in the endoplasmic reticulum were highly enriched. This result indicated that N limitation might have a significant effect on the methylation of carbon-N-metabolism-associated genes. In addition, the pathways involving transcription-related pyrimidine and RNA polymerase were also overaccumulated. In the roots, the terms of protein export, ubiquitin-mediated proteolysis, and proteasome and amino acid (glutamate, etc.) metabolism were largely enriched. It suggested that DNA methylation might occur on the genes involved in organic N-derived N recycling. Therefore, we focused on the effect of N limitation on DNA methylation of some N-metabolism associated genes. Among the genome-wide NRT1.1/NPF6.3 homologs and nine high-affinity NRT2.1 homologs in the rapeseed genome [28], we found only a CHH hypomethylation site located in the promoter of an NRT1.1 homolog (BnaC08g15370D) and a CHH hypomethylation site located in the promoter region of an NRT2.1 homolog (BnaC06g04560D) in the roots under N limitation, respectively (Table 1). N starvation induced an increase in the expression of all the NRT1.1/NPF6.3 and NRT2.1 homologs ( Figure 2E), which might be correlated with their hypomethylation in their promoter regions. However, we also found that DNA methylation in promoter regions did not consistently repress gene expression levels. For instance, we identified a hypermethylated site in the promoter of NRT1.5/NPF7.3 (BnaC05g24580D) in the roots and a hypermethylated site in the promoter region of NRT1.7/NPF2.13 (BnaC06g30920D) in the shoots of rapeseed plants under N limitation (Table 1), whereas their expression was induced by N limitation ( Figure 2E). Furthermore, some other NPF family genes, such as NRT1.6/NPF2.12, NRT1.9/NPF2.9, NRT1.11/NPF1.2 and NAXT1 (nitrate excretion transporter 1)/NPF2.7, also showed significantly differential methylation between N sufficiency and N limitation (Table 1). CLCa, a chloride channel gene, is mainly responsible for vacuolar nitrate influx [29], and we found a hypermethylated site in its promoter region (Table 1), which was consistent with its reduced expression under N limitation ( Figure 2E).
In addition to the molecular identification of NRTs, we also tested a hypermethylated site in an intron of Nitrate Reductase 1 (NIA1, BnaA02g18610D) and a hypermethylated site in an exon of Glutamate Synthetase 1;1 (Gln1;1, BnaC04g30250D) in the roots under N limitation (Table 1), which might be associated with their enhanced transcriptional expression [22]. These above findings indicated that N limitation had a significant effect on N-metabolism-associated genes, including those responsible for N absorption, translocation, remobilization, and assimilation.

Transcriptional Regulation of Rapeseed miRNAs in Response to N Limitation
To determine molecular responses of miRNAs to N limitation and identify the genes targeted by differentially expressed miRNAs, we constructed a total of 12 small RNA libraries from N-sufficient shoots (shoot_0 h), N-sufficient roots (root_0 h), N-deficient shoots (shoot_72 h), and N-deficient roots (root_72 h) for high-throughput sequencing of small RNAs. After removing low-quality reads and adapter sequences, we obtained about 16 million clean reads, corresponding to~0.87 Gb of sequencing data, for each sample (Supplementary Table S2). Based on the results of normalized expression between each two biological replicates, the Pearson correlation coefficients (R) were calculated, most of which were more than 0.92 (Supplementary Figure S8). Low base calling rates, similar GC content, and high correlation among different biological replicates indicated high quality of the sequencing data.
Length distribution of miRNAs frequently reflects the specificity of a specific species or tissue [30]. We found that the nucleotide (nt) lengths of Bna-miRNAs ranged from 18 nts to 30 nts, and the miRNAs with 24 nts were most abundant in all the miRNA libraries ( Figure 6A). Distribution characteristics of the miRNA nucleotide lengths were similar to those of previous reports in B. napus [31,32]. Most Bna-miRNA sequences showed a strong bias for a uridine (U) at the first nucleotide position (Supplementary Figure S9A, Supplementary Table S3). Before the N-limitation treatment, namely at 0 h (N sufficiency), most of the 24-nt long miRNAs had an apparent preference for U in both shoots and roots; however, at 72 h (N limitation), they had a bias for guanine (G) (Supplementary Figure S9B). The hierarchical clustering and heat maps of genome-wide miRNA expression revealed that similar expression patterns occurred in the three biological replicates of each sample. However, significantly differential expression profiling of miRNAs was identified among different plant tissues or N treatments (Supplementary Figure S10).  In this study, the miRNAs showing high expression levels, accounting for a major proportion of the genome-wide miRNAs of rapeseed (Supplementary Table S4), were easily detected in the subsequent studies. We found that different members in the same miRNA family exhibited similar expression levels, such as miR164s and miR167s in the shoots and miR156s and miR171s in the roots. A total of 123 and 101 differentially expressed miRNAs in the shoots and roots, respectively, were identified, more than half of which were downregulated by low N ( Figure 6B). Furthermore, through a Venn diagram analysis, we found that 20 and 42 differentially expressed miRNAs were upregulated and downregulated in both shoots and roots, respectively, ( Figure 6C). To identify biological function of the genes targeted by differentially expressed miRNAs in B. napus during exposure to N deficiency, we plotted GO enrichment and KEGG pathways in which they were involved. We found that most potential targets of the conserved miRNAs were transcription factor genes, such as MYBs and Auxin-responsive factors (ARFs). For example, both miR160 and miR167 target ARF transcription factors [33,34], and miR156 and miR162 target squamosal promoter-binding-like (SPL) proteins and Dicer-liker (DCL) proteins. In addition, target genes that are implicated in leaf development and cell division, as well as phytohormone signaling, were also highly accumulated ( Figure 6D). In our previous study, we identified numerous MYB transcription factor genes that were differentially expressed between N sufficiency and N limitation [22], which might be mainly involved in the anthocyanin biosynthesis in response to N limitation. For the KEGG analysis, the pathways involving environmental adaptation, carbon (such as carbohydrate and lipid) and N (such as amino acid) metabolism, and membrane transport were significantly accumulated ( Figure 6E).
In general, miRNAs are highly conserved among different plant species [35]. We used the miRbase database [36] to detect known miRNAs that were differentially responsive to high N and low N conditions. A total of 34 known differentially expressed miRNAs belonging to 22 families and 89 novel were identified in the shoots ( Figure 7A,B). We found a total of 22 known differentially expressed miRNAs belonging to 16 families and 79 novel differentially expressed miRNAs in the roots ( Figure 7C,D). Among the known miRNAs, the Bna-miR169/Bna-miR167 families in the shoots and the Bna-miR156/Bna-miR171 families in the roots had the most miRNA members, and the Bna-miR171 family members presented the largest differential expression between N sufficiency and N limitation in both shoots and roots ( Figure 7A,C). In the shoots, both Bna-miR164a and Bna-miR164b, which regulated leaf senescence, were significantly induced under N limitation ( Figure 7A). Bna-miR156, which regulated production of anthocyanin, was downregulated in the shoots and was upregulated in the roots [37]. In the roots, both Bna-miR160a and Bna-miR171 family members (including miR171a, miR171f, and miR171g), which regulated root growth [38], were significantly elevated under N limitation compared with those under N sufficiency ( Figure 7C), and they might facilitate N-starvation-induced root elongation. Overexpression of miR169 repressed the expression of the major NRT genes AtNRT1.1 and AtNRT2.1 [15].
In this study, the expression of Bna-mi169m was decreased in the roots under N limitation ( Figure 7C), which might be associated with the enhanced transcript levels of BnaNRT1.1s and BnaNRT2.1s ( Figure 2E). Furthermore, downregulation of miR164 combined with the upregulation of NAC1 produced more lateral roots [39]. However, we found that two miR164 homologs (Bna-miR164a/b) were upregulated in the shoots and roots under N limitation ( Figure 7A,C). Therefore, we presumed that Bna-miR164a/b was not involved in the N-starvation-induced lateral root proliferation.

Degradome Sequencing-Assisted Identification of Genes Targeted by Differentially Expressed miRNAs in Response to N Limitation
To identify miRNA-directed mRNA cleavage events in rapeseed plants under N limitation, we subjected rapeseed plants under high N and low N to a high-throughput degradome sequencing platform. A total of 18,370,535 clean reads of 23.67 million raw reads in the shoots and 20,131,661 clean reads of 20.23 million raw reads in the roots were obtained, and more than 70% of the clean reads were matched with the reference sequences assembled from the rapeseed transcriptome database (Supplementary Table S5). In total, we identified 31 genes identified as targets of known miRNAs, which covered 30 miRNA families, consisting of 48 members (Supplementary Table S6).
Under N limitation, we found not only chlorophyll decrease, anthocyanin accumulation, leaf branch reduction, and leaf senescence in the shoots but also induction of primary and lateral roots (Figure 1), which were considered as adaptive nutrient morphogenesis responses of rapeseed plants to N limitation [22]. To identify the miRNA-mediated regulatory mechanisms underlying these morpho-physiological responses, we tested the expression of several related miRNAs and their target genes. In our present and previous studies, we identified increased expression of Bna-miR827 in the shoots ( Figure 7A). The involvement of the At-miR827-AtNLA-AtORE1 cascade pathway in the N-limitation-induced leaf senescence in Arabidopsis was validated [17]. We found a higher transcript accumulation of Bna-miR827 in the leaves of rapeseed plants under N limitation than that under N sufficiency ( Figure 7A), which repressed the expression of NLA [22] and further led to the enhanced expression of ORE1, a major transcription factor regulating leaf senescence.
In Arabidopsis, SCL6 (SCARECROW-Like 6), which is targeted by miR171, plays an important role in the regulation of shoot branch production and leaf morphogenesis [40]. In this study, we found that the identified Bna-miR171f was induced in both shoots and roots under N limitation ( Figure 7A,C). The mature sequence of Bna-miR171f was mapped onto the 3'-end of its precursor ( Figure 8A), and it was highly identical with the corresponding homolog in Arabidopsis ( Figure 8B). Stem-loop RT-qPCR confirmed that the expression of Bna-miR171f was induced by N limitation in both shoots and roots ( Figure 8C). Furthermore, we examined the expression of the SCL6 homologs in B. napus and found that only BnaAn.SCL6b (BnaAnn18550D) was significantly repressed by N limitation ( Figure 8D). The antagonistic expression between Bna-miR171f and BnaAn.SCL6b indicated their potential interaction. Subsequently, we used degradome sequencing to further validate the presumption. The degradome result showed that the Bna-miR171f-mediated cleavage site fully matched the transcript of BnaAn.SCL6b, whose transcript had only an obvious peak showing the cleavage site ( Figure 8E,F).
In A. thaliana, overexpression of ARF17, which is targeted by miR160, leads to root length reduction [41]. In this study, we identified that Bna-miR160a, induced by N limitation (Figure 7A), was mapped onto the 5'-end of its precursor ( Figure 9A), and that it shared a high identity with its homolog in Arabidopsis ( Figure 9B). RT-qPCR confirmed that the expression of Bna-miR160a was induced by N limitation in both shoots and roots ( Figure 9C). Among the ARF17 homologs in B. napus, only BnaC6.ARF17a (BnaC06g20640D) was differentially expressed between high N and low N conditions ( Figure 9D). Furthermore, the degradome data showed that the Bna-miR160a-mediated cleavage site fully matched the transcript of BnaC6.ARF17a, whose transcript had only an obvious peak showing the cleavage site ( Figure 9E,F).

Discussion
N nutrients play pivotal roles in modulating crop growth, organ development, and seed productivity [42]. Currently, great process has been made in uncovering the genetic and molecular mechanisms underlying NUE and NLA in the model plants Arabidopsis and rice [43][44][45][46][47][48][49]. In addition to characterization of the aforementioned functional structural genes, epigenetic modification, such as phosphorylation of NRT1.1 [50], methylation of NRT2.1 [51], and post-transcriptional cleavage of NLA by miR827 [17], is also significant for the regulation of NUE and NLA in plants.
Rapeseed has a higher demand for N nutrients than do many other crop species [20,52]. For the purpose of guaranteeing seed yield production, enhancing rapeseed NLA is a core route for not only reducing N fertilizer applications but also increasing NUE. However, the implicated mechanisms underlying NUE and NLA in allotetraploid rapeseed remain unclear. In our study and other previous studies, genome-wide transcriptomic and proteomic responses of rapeseed plants to N limitation were investigated through high-throughput sequencing [5,22]. Furthermore, in this study, we employed WGBS and small RNA sequencing to identify epigenetic modification, including DNA methylation and miRNA cleavage, in regulating adaptive responses of rapeseed plants to N limitation.

N Limitation Triggered Adaptation and NUE Enhancement
Nitrate functions not only as a N nutrient source, but as a signal molecule affecting plant growth and development in plants by inducing expression of auxin-related genes [53], modulating leaf growth and development [54,55], and regulating root system architecture [56][57][58][59][60]. An Arabidopsis mutant, nla, shows leaf early senescence and does not accumulate anthocyanin under N limitation, not exhibiting an adaptive response of plants to limited N nutrients [17]. Under N limitation, the rapeseed leaves and hypocotyls rapidly accumulated large quantities of anthocyanins ( Figure 1A). Anthocyanins belong to important flavonoid categories, and they can protect plants from N-deficiency-induced damages of reactive oxygen species [61]. Under N limitation, the NUE of rapeseed plants was significantly increased ( Figure 2C). Startup of limited N induces leaf senescence, which triggers N nutrient remobilization from senescent leaves to developing leaves or other organs, such as seeds [17]. However, weak NLA induces leaf early senescence, easily resulting in leaf dropping and detaching from plants before that N can be sufficiently remobilized to the sink organs [17]. Therefore, it suggested that improving plant adaptation to N limitation is an urgent route for NUE enhancement.
As a major determinant of nutrient uptake, developmental plasticity of RSA has a significant impact on plant growth performance and nutrient use efficiency in nutrient-depleted soils [62]. RSA-based genetic improvement is proposed to be a feasible tool for the rapeseed NUE enhancement [63]. In turn, RSA is also influenced by N availability, which allows plants to optimize N acquisition from growth medium by modulating RSA [64]. In this study, we found that under N deficiency (0.30 mM NO 3 − ), the primary roots of rapeseed became longer, combined with more lateral roots than those under N sufficiency (6.0 mM NO 3 − ) ( Figure 1H). It might contribute to root foraging for N nutrients in the growth media ( Figure 10). Based on these findings, we proposed that identification and characterization of the genes regulating RSA is urgent for the genetic improvement of plant NLA and NUE.

Genome-Wide Differential DNA Methylation Reveals Pivotal Roles of Epigenetic Regulation in Adaptive Responses of Rapeseed to N Limitation
DNA methylation plays an important role in regulation of transposon silencing, genome stability, and gene expression in plants, and its pattern is affected by environmental biotic and abiotic factors. Nutrient stresses in soils, such as phosphate starvation and N deficiency, were shown to alter DNA methylation status within plants [8,65]. Prior to this study, few research studies were conducted to identify and characterize the DNA methylation fingerprints of allotetraploid rapeseed under N limitation. Therefore, to illustrate the effect of N limitation on DNA methylation, we showed how rapeseed plants remodeled methylation landscapes of genomic DNA across gene structures in response to N limitation. , all of which contributed to stronger root system that facilitated efficient N forage. Limited N supply repressed the expression of CLCa through promoter hypermethylation regulation, and it reduced vacuolar N storage, which was favorable for N assimilation in the cytoplasm. In addition, N limitation also enhanced the NRT1.5 expression and repressed the NRT1.8 expression, both of which facilitated N translocation from roots to shoots. Furthermore, in the shoots, limitation of N induced an increase in the miR827 expression that downregulated the E3 ubiquitin ligase NLA, which accelerated the remobilization of N nutrients from old leaves to young ones through modulating NRT1.7 degradation. The hypomethylation of the nitrate reductase NIA1 and hypermethylation of Gln1;1 within their gene bodies might be responsible for the reduced expression of NIA1 and increased expression of Gln1;1. Green solid arrows and yellow solid lines indicate the valid processes or genes are induced and repressed by N limitation, and the dashed arrows indicate potential DNA methylation under N limitation. Red words indicate the epigenetic regulation of N-metabolism genes identified in this study, and grey words refer to regulatory mechanisms that were previously reported.
Despite derivation from a common ancestor, the diploid Brassica oleracea (C o C o ) and Brassica rapa (A r A r ) contain 38.80% and 21.47% transposable elements (TEs), respectively, in their assembled genomes, respectively [66]. The above facts lead to an asymmetric distribution of TEs in the A n and C n subgenomes of the assembled B. napus genome (A n A n C n C n ), which may lead to different methylation patterns between the two subgenomes [19]. Regulatory roles of DNA methylation in gene transcription are dependent on plant tissue types, methylated genomic contexts, and target gene regions [67]. DNA methylation of regulatory elements, such as gene promoters and enhancers, potentially repress gene expression by altering chromatin structures and inhibiting transcription initiation [68]. In this study, we identified hypomethylation in the gene promoter regions of rapeseed NRT1.1 and NRT2.1 (Table 1), which might be involved in their transcription activation under N limitation ( Figure 2E). A hypermethylated site in an exon of Gln1;1/GS1;1 (BnaC04g30250D) was identified in the roots of rapeseed plants under N limitation (Table 1), whereas its expression was enhanced ( Figure 2E). Therefore, gene body methylation does not necessarily reduce the gene expression [69], but instead plays a positive role in transcription activation in some plant species, such as Arabidopsis [70], rice, and maize [71]. Furthermore, we found that DNA methylation patterns of some genes were not consistent with their transcriptional expression, such as NRT1.5 and NRT1.7 ( Figure 2E, Table 1). MYB59 directly bound to the NRT1.5 promoter and positively regulated its expression [72]. NLA which encodes a RING E3 ubiquitin ligase, interacts directly with NRT1.7 in the plasma membrane and promotes NRT1.7 degradation [17]. The above results indicated that other regulatory mechanisms-for instance, transcription factor binding or post-translational modification-might be involved in the transcriptional modulation of these genes above-mentioned.

Global miRNA and Degradome Sequencing-Assisted Identification of Differentially Expressed miRNAs Uncover Epigenetic Regulatory Mechanisms Underlying N-Limitation-Induced Morphogenesis
To identify the roles of mRNA, miRNA, and their targets in regulating these N-deficiency-induced morphogenesis, we employed high-throughput genome-wide mRNA sequencing, miRNA sequencing, degradome sequencing, and RT-qPCR assays to identify differential expression profiling in response to N limitation. Genome-wide mRNA transcriptome sequencing and RT-qPCR assays revealed that under N limitation, the expression of NRT1.1, NRT1.5, NRT1.7, NRT2.1/NAR2.1, and Gln1;1 was significantly increased, whereas the transcript levels of CLCa, NRT1.8, and NIA1 were reduced ( Figure 2E). This finding indicated that limited N improved plant NUE by enhancing efficient N absorption, reducing vacuolar N storage and root unloading, facilitating root-to-shoot long-distance translocation, improving source-to-sink remobilization, and strengthening N assimilation ( Figure 10).
Auxin-mediated nitrate signaling by NRT1.1 participates in adaptive response of Arabidopsis RSA to the spatial heterogeneity of nitrate availability [73]. Moreover, NRT2.1 is also shown to be involved in nitrate-dependent root elongation by regulating auxin transport to roots [74]. Enhanced expression of nitrate transporters plays a key role in the regulation of RSA-mediated adaptive responses to N limitation through modulation of auxin signaling. Furthermore, we employed miRNA and degradome sequencing to identify key miRNAs and their targets that were involved in regulating adaptive responses of rapeseed plants to N limitation. In general, we found that the genes that are involved in leaf development and phytohormone (particularly auxin) signaling were enriched ( Figure 6D,E). Therefore, we presumed that differentially expressed miRNAs might be implicated in the regulation of N-limitation-induced shoot and root morphogenesis, and we proposed that a genetic engineering approach by manipulation of miRNA expression may be a feasible route for the improvement of NUE and NLA in rapeseed plants. Uniform rapeseed seedlings were cultivated in an illuminated growth chamber (300-320 µmol m −2 s −1 ; 24 • C daytime/22 • C night; 16 h light/8 h dark) with relative humidity of 75%. For nitrate-depletion treatments, rapeseed seedlings after 5-day germination were first hydroponically grown under high nitrate (6.0 mM) for 10 days, and then the rapeseed plants were transferred to low nitrate (0.30 mM). At the times of 0 h and 72 h after N limitation treatment, rapeseed shoots and roots were individually sampled for isolation of genomic DNA and total RNA. The solution concentration of N nutrients was adjusted by reducing KNO 3 and replacing Ca(NO 3 ) 2 with CaCl 2 , and K + was complemented by adding KCl [63].

Quantification of Morpho-Physiological Characteristics
Total leaf areas of rapeseed plants were determined using an LI-COR LI-3100C leaf area meter. Anthocyanin was isolated using extraction buffer (CH 3 OH:H 2 O:H 2 O 2 = 60:13:2, v/v/v); then, the purified extract was assayed by a UV-1800 spectrophotometer on the absorbance wavelengths at 530 and 657 nm [22]. The roots of rapeseed seedlings cultivated under hydroponic culture system were subjected to a image scanner, and then total root length, root volume, root surface area, and root average diameter were analyzed through WinRHIZO Pro (Regent Instruments, QC, Canada) [75]. Soluble sugars were isolated as described by [21], and were assayed for glucose, fructose, and sucrose concentrations using a commercially available kit (Megazyme, http://www.megazyme.com/).
Harvested rapeseed plants were oven-dried at 65 • C until a constant weight was obtained, and then the dried samples were digested in a HNO 3 /HClO 4 (4:1, v/v) mixture at 200 • C until completely digested. Concentrations of mineral elements were quantified by ICP-MS (NexION TM 350X; PerkinElmer, Massachusetts, USA). Fresh tissues of rapeseed plants were sampled and bathed in boiling water for 30 min; then, they were assayed for nitrate concentrations using a continuous flow autoanalyzer (AA3, Seal Analytical Inc., Southampton, UK) [76]. After treatment at 120 • C for 30 min, the samples were dried at 65 • C for 72 h, and then were milled to a fine powder. Total N concentrations of oven-dried plant tissues were determined with a carbon and N analyzer (Elemental Analyzer EA1108; Carlo Erba Strumentazione, Milan, Italy). NUE values were calculated as the following formula: NUE = total plant dry weight (g)/total plant N accumulation (g) [77].

DNA and RNA Isolation, Quantification, and Qualification
Genomic DNA (gDNA) and total RNA was isolated using 2% CTAB and TRIzol, respectively. A Nanophotometer ® spectrophotometer (IMPLEN, CA, USA) was used to check the purity of DNA and RNA, and their concentrations were measured using a Qubit ® DNA and RNA Assay Kit in a Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA). An RNA Nano 6000 Assay Kit of an Agilent Bioanalyser 2000 system (Agilent Technologies, CA, USA) was used to assess RNA integrity by determining the RNA integrity number.

Library Preparation for Whole Genome Bisulfite Sequencing (WGBS) and miRNA Sequencing
A total amount of 5.0 µg gDNA spiked with λDNA was fragmented by sonication to 200-300 bp with a Covaris S220, which followed by end repair and adenylation. Cytosine-methylated barcodes were ligated to the sonicated DNA. The DNA fragments were treated twice with bisulfite using EZ DNA Methylation-GoldTM Kit (Zymo Research) before the ssDNA fragments were amplified by PCR using KAPA HiFi HotStart Uracil + ReadyMix (2X). Library concentration was quantified by Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA) and quantitative PCR, and the insert size was assayed on an Agilent Bioanalyser 2100 system (Agilent technology).
A total of about 3.0 µg RNA per sample was used for construction of small RNA libraries, which was generated using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB, USA). Briefly, the NEB small RNA adapters were ligated to the 3 and 5 ends of miRNAs. Reverse transcription primers were hybridized to the excess of 3 small RNA adapter, and the first strand DNA was synthesized using the M-MuLV reverse transcriptase (RNase H − ). PCR amplification was performed using LongAmp Taq 2X Master Mix, small RNA primers from Illumina, and index primers. PCR products were purified through 8% polyacrylamide gel electrophoresis, and the DNA fragments corresponding to 140-160 bp were recovered. Library quality was assessed on an Agilent Bioanalyser 2100 system (Agilent Technologies, CA, USA).

Clustering and Sequencing
For WGBS, the DNA libraries were sequenced on an Illumina Hiseq 4000, and 150 bp paired-end reads were generated.
For the miRNA sequencing, the clustering of index-code samples was performed on a cBot Cluster Generation system using TruSeq SR Kit v3-cBOt-HS (Illumina). After cluster generation, the libraries were sequenced on an Illumina Hiseq 2500 platform, which generated 50-bp single-end reads.

Quality Control
Read sequences (FASTq) produced by the Illumina pipeline were preprocessed through Trimmomatic (v. 0.35) as follows: first, reads containing adapter sequences were filtered out; second, some reads having ambiguous bases in their sequences were removed; third, reads with more than 50% low quality bases (PHRED score ≤ 20) were excluded. Subsequently, Q 20 , Q 30 , and GC content of clean data were calculated for further analysis.

Reads Mapping to the Reference Sequences
Bismark (v. 0.12.5) was used to perform alignments of bisulfite-treated reads to the reference genome of B. napus (ZS11) annotated in the BnPIR (Brassica napus pan-genome information resource) (http://cbi.hzau.edu.cn/bnapus/index.php) database. The reference genome was first transformed into the bisulfite-converted version, and then were indexed using Bowtie2. Sequence reads were also transformed into the fully bisulfite-converted versions (C-to-T and G-to-A conversion) before they were aligned to the similarly converted versions of the genome in a directional manner. Sequence reads that produced a unique best alignment from the two alignment processes (original top and bottom strand) were then compared with the normal genomic sequences, and the methylation states of all cytosine positions in the read were inferred. Read pairs that shared the same coordinates in the rapeseed genome were regarded as duplicated and were removed before methylation state determination to avoid potential calculation bias of the methylation level. The sodium bisulfite non-conversion rate was calculated as the percentage of cytosine sequenced at cytosine reference positions in the lambda genome.

miRNA Alignment and Target Genes Prediction
The miRBase 22.1 (http://www.mirbase.org/) [79] database was used to retrieve published miRNA sequences and functional annotation. miRDeep2 was used for the identification of novel and known miRNAs. The miREvo and miRDeep2 databases were integrated to predict novel miRNAs through exploring secondary structures, Dicer cleavage sites, and minimum free energy of small RNA tags without annotation. Prediction of target genes of miRNAs was performed by psRobot_tar in psRobot.

Differential Analysis of DNA Methylation and miRNA Expression
After removing low-quality reads (Q score < 28), 150 bp paired-end WGBS reads were mapped to the rapeseed genome using Bismark software. Only cytosines that were covered by at least three reads were considered and counted [80]. The methylation level of individual cytosines was calculated as the ratio of mC to total cytosines [mC/(mC+un-mC)]. Fold change at each site was calculated by comparing the methylation ratio under N deficiency to the methylation ratio under N sufficiency. An mC fold change greater than 1.2 was considered to indicate differential methylation. Densities and average methylation levels were calculated in 50 kb units. The Kruskal-Wallis rank sum test was performed on each bin for the mCs in CG, CHG, and CHH contexts to identify significant differences between N sufficiency and N deficiency. The p-value of the bins was adjusted using BH multiple test correction. Bins with corrected p < 0.05 were defined as differentially methylated regions (DMRs). Differentially methylated genes (DMGs) were defined as genes whose gene body or promoter region (upstream 2.0 kb from the transcription starting site) have at least 1 bp overlap with their DMRs. Hypermethylation/hypomethylation status was assigned to each DMR based on the differences in methylation levels between the two N levels. If the methylation level of a DMR under N sufficiency was lower than that under N limitation, the DMR was considered a hypermethylated site and vice versa [81]. miRNA expression levels were estimated by transcripts per million (TPM). Corrected p ≤ 0.05 and |log 2 (fold change)| ≥ 0.60 were set as the threshold for the identification of differentially expressed miRNAs using the DESeq R package (1.10.1).

Gene Ontology and KEGG Enrichment Analysis
Differentially expressed genes (cutoff: p < 0.05 and false discovery rate < 0.05) that were identified by transcriptome sequencing were submitted to the online Protein ANalysis THrough Evolutionary Relationships (PANTHER) classification system (http://www.pantherdb.org/) for Gene Ontology (GO) enrichment analysis [82].
In terms of the data of WGBS and miRNA sequencing, the GOseq R package based on Wallenius non-central hyper-geometric distribution was used to perform GO enrichment analysis of the target genes. DEG enrichment in Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.kegg.jp/) pathways was assessed using the KOBAS 3.0 software (http://kobas.cbi.pku.edu.cn/kobas3) [83,84]. Both p-values and false discovery rate values less than 0.05 were used to identify overrepresentation of GO categories and KEGG pathways.

Degradome Analysis
For degradome sequencing, the rapeseed seedlings, ZS11, which were first hydroponically grown under high nitrate (6.0 mM) for 7 days, were then transferred to low nitrate (0.30 mM). At the times of 0 h and 72 h after N limitation treatment, rapeseed shoots and roots were individually sampled for RNA isolation. Parallel analysis of RNA ends libraries was performed on an Illumina's Cluster Station and sequencing was subsequently conducted on Illumina GAIIx. CLEAVELAND 3.0 was used for the sequencing data analysis.

Reverse Transcription Quantitative PCR (RT-qPCR) Assays
Shoots and roots of rapeseed plants were individually harvested and immediately stored at −80 • C until RNA isolation. Each sample contained three independent biological replicates for transcriptional analyses. Reverse transcription quantitative PCR (RT-qPCR) assays were used to determine the gene expression. After removing genomic DNA from the RNA samples with RNase-free DNase I, complementary DNA (cDNA) synthesis was performed using the PrimeScript TM RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Shiga, Japan) with total RNA as the templates. Quantitative analysis of relative gene expression was performed using the SYBR ® Premix Ex Taq TM II (Tli RNaseH Plus) (TaKaRa, Shiga, Japan) kit in an Applied Biosystems StepOne TM Plus Real-time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). Thermal recycle regimes were as follows: 95 • C for 3 min, followed by 40 cycles of 95 • C for 10 s, then 60 • C for 30 s. A melting curve analysis was also conducted to ensure primer specificity of target genes: 95 • C for 15 s, 60 • C for 1 min, and 60-95 • C for 15 s (+0.3 • C per cycle) [85]. Expression data were normalized using the data of public reference genes BnaEF1-α [86] and BnaGDI1 [87], and relative gene expression was calculated with the 2 −∆∆C T method [88]. Stem-loop reverse transcription PCR was used in the RT-qPCR experiments for miRNAs [89], and BnaU6 snRNA was used as an internal control in each reaction [90]. Sequences and amplification efficiency of primers used in this study are listed in additional file: Supplementary  Table S7.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1422-0067/21/22/8453/s1, Figure S1: Overview of the genome-wide differential DNA methylation fingerprints in the roots of rapeseed plants    The funding organizations provided the financial support to the research projects, but were not involved in the design of the study, data collection, analysis of the data, or the writing of the manuscript.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.