MiRNA-mediated Changes in DNA Methylation Affect the Expression of Genes Involved in the Thickness of Pod Canopy Trait in Brassica Napus

Background: Methylation plays an important role in regulating crop development, but little is known about how methylation regulates plant architecture in rapeseed (Brassica napus). Here, we examined how methylation affects the TPC (thickness of pod canopy) trait in rapeseed by performing genome-wide methylation analysis of two extreme TPC lines. Results: We detected signicant differences in overall methylation levels between the high- and low-TPC lines in the CG, CHG, and CHH contexts in the promoters of genes in the stem apex and ower bud. In ower buds, 26 genes had signicantly higher methylation levels in the high-TPC samples compared to the low-TPC samples, resulting in signicantly reduced gene expression. By contrast, in the stem apex samples, the promoter regions of 22 genes were hypermethylated in the high- vs. low-TPC samples. The promoters of 19 and 21 genes had signicantly reduced methylation levels in the ower bud and stem apex, respectively, of the high- vs. low-TPC samples, resulting in signicantly higher expression levels. Some of these differentially expressed genes are associated with TPC-related traits, such as BnaC03g53050D (UBC32), BnaA05g26660D (CYSB), BnaA10g07880D (TCP 1), BnaAnng09670D (SMP1), BnaA09g02000D (SDH2-2), BnaC01g12960D (NRT1.8), and BnaC09g30490D (TAF15b). In addition, 14 important genes related to growth and development were differentially regulated between the two groups due to miRNA-mediated differences in methylation levels in their promoters. For example, hypermethylation in the promoter region of BnaCnng64040D (Lipase family protein) mediated by miR159a led to signicantly reduced gene expression in ower buds of high-TPC vs. low-TPC lines. Conclusions: These results, together with our previously generated RNA-seq and miRNA proling data, indicate that both methylation and miRNAs are involved in regulating the expression of genes in nitrogen-related metabolic pathways, thereby affecting the TPC trait in B. napus, providing a reference for uncovering the molecular mechanism regulating this crucial trait. the expression of BnaCnng64040D (Lipase family protein) regulated by miR159a, hypermethylation in its promoter region to a signicant decrease in gene expression in the ower buds of high- vs. low-TPC samples. The expression of BnaC09g30490D (TAF15b) and BnaC03g09180D is regulated by miR167c and miR827, respectively, and the loss of methylation in their promoter regions resulted in a signicant increase in their expression in the ower buds of the high-TPC lines. Both BnaC02g22120D (Spc97/Spc98 family of spindle pole body component) and BnaA07g07800D are regulated by miR319. Increased promoter methylation led to a signicant decrease in the expression of these genes in the stem apex of the high-TPC lines. BnaA04g14930D (S-locus lectin protein kinase family protein) and BnaCnng56050D (Pyridoxal phosphate-dependent Transferases superfamily protein) are regulated by miR159, miR319, miR122, and so on. The loss of methylation in the promoters of these genes led to a signicant increase in gene expression in the stem apex and ower buds of the high-TPC lines.


Introduction
Oilseed rape (Brassica napus, 2n = 38, AACC) is one of the most important oilseed crops worldwide. Yield improvement is an extremely important goal for all rapeseed breeders. Modifying plant type is the most important way to increase rapeseed yields. Thickness of pod canopy (TPC), i.e., the thickness of the plant canopy from the bottom-most effective (seed-bearing) pod to the uppermost effective pod, is a key trait that determines the three-dimensional structure of plants. TPC is directly related to other important traits, such as economic yield (EY, grain weight per plant), plant height (PH), pod terminal height (PTH), stem height (SH), rst effective branch height (FEBH), height of the lowest pod (HLP), rst effective branch number (FEBN), main in orescence effective length (MIEL), and rst uneffective branch number (FUBN) [1]. Hence, understanding the molecular mechanism of TPC could shed light on the formation of plant architecture and help breeders further increase production in B. napus.
In recent years, many studies have focused on traits in crops. PH is an important factor affecting rice yield. OsMPH1 improves grain yield by regulating PH in rice [2]. MiR319 expression induces dwar sm to suppress PH in rice [3]. A genome-wide association study (GWAS) of PH and primary branch number in rapeseed revealed eight PH-related quantitative trait loci (QTLs) on chromosomes A03, A05, A07, and C07 and ve PBrelated QTLs on chromosomes A01, A03, A07, and C07 [4]. GWAS identi ed four, ve, and seven SNPs for FEBH, FEBN, and PH, respectively, in B. napus and uncovered many genes associated with these traits [5]. We recently reported that TPC is regulated by miRNAs and identi ed many candidate genes for TPC based on RNA-Seq and miRNA pro ling analyses [1]. Many genes involved in nitrogen-related responses were dramatically differentially expressed in high-vs. low-TPC lines, such as ASP5, ASP2, ASN3, ATCYSC1, PAL2, APT2, CRTISO, and COX15 [1], suggesting that nitrogen metabolism plays important roles in regulating TPC in B. napus.
DNA methylation is a critical epigenetic modi cation in many plants. The methylation and demethylation of genes are dynamically regulated during plant growth [6]. Plant DNA is methylated in three sequence contexts: symmetrical mCG and mCHG and asymmetrical mCHH (H = A, T, or C) [7][8].
Methylation is regulated and maintained by independent pathways that can be broadly classi ed into maintenance of methylation and de novo methylation pathways [6]. Due to their symmetrical nature, CG and CHG sites are substrates for maintenance methyltransferases, which recognize hemimethylated DNA after replication and methylate newly synthesized unmethylated strands [9]. De novo methyltransferases establish DNA methylation across all sequence contexts and are required for maintaining asymmetric (CHH) methylation.
Small RNAs and proteins related to modi ed histones affect DNA methylation [6]. DNA methylation and demethylation are distinct process: the latter is regulated by 5-methylcytosine DNA glycosylase enzymes in plants [10]. 5-methylcytosine is removed from DNA by base excision repair, which is essential for the expression of imprinted genes and endosperm development [11].
In plants, DNA methylation plays a key role in the activation of maturation-inducing genes and the inhibition of maturation-inhibiting genes in tomato (Solanum lycopersicum) [12]. DNA methylation increases substantially during citrus (Citrus sp.) fruit development and maturation [13]. DNA methylation is closely related to the development of many grain crops. Genome-wide changes in DNA methylation are related to drought stress tolerance [14] and cadmium stress in plants [15]. Overall DNA methylation levels increase during somatic embryogenesis in soybean (Glycine max) [16]. Cytosine methylation plays a positive role in regulating iso avone synthase gene expression and iso avonoid biosynthesis in soybean seeds [17]. Low zinc levels lead to the loss of methylation in maize (Zea mays) roots; conversely, the loss of methylation in root cells leads to zinc de ciency [18]. DNA methylation changes dynamically in response to heat stress in maize seedlings [19]. Lead, cadmium, and zinc toxicity alter DNA methylation levels, thereby enhancing heavy metal tolerance in wheat (Triticum aestivum) [20]. Genome-wide DNA methylation pro ling of ower buds revealed the role of DNA methylation in the molecular regulation of genic male sterility in B. napus [21]. Finally, short-term heat-shock treatment of cultured B. napus microspores led to global changes in DNA methylation, indicating that DNA methylation plays a key role in heat-stress responses [22].
Methylation also affects plant architecture. The proper H3K4me3 levels, which are regulated by COMPASS-like complexes, are critical for rice development and can affect owering and branching [23]. In addition, changes in methylation at the rice (Oryza sativa, Os) FERTILIZATION INDEPENDENT ENDOSPERM 2 locus play important roles in regulating plant height and yield in rice [24]. Changes in methylation are also related to nal plant height in Arabidopsis thaliana [25]. However, our understanding of the mechanisms underlying how methylation regulates plant-type traits in crops remains limited.
Small RNA-directed DNA methylation (RdDM) is a key regulatory pathway that affects numerous plant traits by altering the expression of various genes. The ARGONAUTE 4 (AGO4)-dependent RdDM pathway represses the expression of HOMOLOG OF RPW8 4 and alters the response to submergence in Arabidopsis (Arabidopsis thaliana) [26]. RdDM inhibits the expression of APETALA3 (AP3); transgenic Arabidopsis plants harboring promoter lines showed abnormal stamens and petals [27]. The RdDM pathway is also involved in regulating the activity of the poly(A) polymerase PAPS1 in Arabidopsis, thereby affecting sporophyte and pollen development [28]. RdDM also regulates seed dormancy in plants [29]. This pathway regulates the expression of numerous genes in crops and affects the corresponding traits. ITRAQ (isobaric tags for relative and absolute quantitation) analysis of the leaf proteome indicated that the RdDM pathway plays an key role in defense against geminivirus-betasatellite infection in tobacco (Nicitiana tabacum) [30]. The reinforcement of DNA methylation at CHH sites regulated by the RdDM pathway leads to decreased gene expression, thereby in uencing somatic embryogenesis in soybean [16]. However, the role of the RdDM pathway in regulating plant-type traits in crops has not been reported.
TPC is an important plant-type trait that is closely associated with various phenotypes. We recently identi ed many genes related to the formation of TPC and identi ed miRNAs that limit their expression in B. napus [1]. In the current study, we performed genome-wide methylation analysis of highand low-TPC B. napus lines. The results of this study, combined with previously obtained RNA-and sRNA-sequencing data, shed light on the molecular mechanism of TPC in this important oilseed crop.

Plant Materials and Bisul te Sequencing
We performed methylation-sequencing analysis of plants cultivated at the Chongqing Rapeseed Engineering Research Center, Southwest University, Chongqing, China (106.40°E, 29.80°N). B. napus lines YC4 (SWU71) and YC33 (10-1047) showed relatively high TPC values based on three years of continuous observation, while YC11 (Zhongshuang11) and YC15 (Zhongyou821) showed relatively low TPC values. We therefore used these lines for comparative analyses. When the plants bolted and produced tiny ower buds, tissue samples from the stem apex and nearby ower buds were collected: this represents the key stage for the formation of TPC (after this stage, the stem begins to elongate and pollen begins to develop). We collected at least ve sets of samples per line, each set from a different individual. Samples were collected at approximately the same time, and all samples from each line were mixed into one pool. Total DNA was isolated from the samples using a Rapid Plant Genomic DNA Isolation Kit (Sangon, China) according to the manufacturer's protocols.
After checking sample quality, we added lambda DNA to the sample as a negative control and used a Covaris S220 sonicator to randomly break the genomic DNA into 200-300 bp fragments. We performed end repair on the broken DNA fragments, added an A tail, and attached them to sequencing adapters whose cytosines were modi ed by methylation. We performed bisul te treatment (using an EZ DNA Methylation Gold Kit, Zymo Research), which converted unmethylated C to U (converted to T after PCR ampli cation), whereas methylated C remained unchanged. PCR was performed to obtain the nal DNA library. A library construction owchart is shown in Supplementary Fig. 1.

Analysis of Sequencing Data
After obtaining the original sequencing reads, we performed bioinformatic analysis as shown in Supplementary Fig. 2. We preprocessed the reads (in FASTQ format) produced by Illumina sequencing using Trimmomatic (version 0.35). The steps were as follows: First, reads containing adapter sequences were ltered out; Second, reads with N (unknown bases) < 10% were deleted; third, reads containing > 50% low-quality bases (PHRED score ≤ 20) were removed. The Q20, Q30, and GC contents of the data were calculated at the same time. The remaining reads that passed all the ltering steps were counted as clean reads and used for subsequent analysis.

Aligning the Reads to the Reference Genome
Bismark software (version 0.12.5) was used to compare the bisul te-treated reads to the reference genome using default values. The reference genome was converted to a bisul te-converted version (C-T and G-A conversion) and indexed by Bowtie2 (version 2.2.5). The reads were also converted into bisul te-converted versions (C-T and G-A conversion) and directly compared with the converted reference genome. The reads were then compared with the normal genomic sequence to infer the methylation status of each cytosine position in the reads. Reads pairs with the same coordinates in the genome were treated as duplicates and deleted before the methylation status was determined to avoid potential calculation bias of the methylation level. The non-conversion rate of bisul te was calculated as the percentage of sequenced cytosine to the reference cytosine in the lambda genome. The basic concept behind Bismark comparison is shown in Supplementary Fig. 3.

Analysis of Differentially Methylated Regions
Differentially methylated regions (DMRs) were identi ed using swDMR software (https://sourceforge.net/projects/swdmr/) with the conditions read coverage ≥ 5, methylation level difference ≥ 0.1 or fold change ≥ 2, corrected p-value < 0.01 using the sliding window method. The window size was set to 1000 bp and the step size was 100 bp. Fisher's exact test was used to detect DMRs. Genes whose functional regions overlapped with DMRs by at least 1 bp were de ned as DMR-associated genes (DMGs).

GO Analysis of DMGs
The DMGs, which were corrected for gene length deviations, were subjected to gene ontology (GO) analysis using the GOseq R package. GO terms with corrected p-values < 0.05 were considered to be signi cantly enriched by the DMGs. All DMGs were annotated with BGI Web Gene Ontology Annotation plot (http://wego.genomics.org.cn/).

FastQC and Trimming of Raw Data
With the rapid development of high-throughput sequencing technology, massive amounts of data can now be obtained by sequencing. Bioinformatic analysis is used to obtain useful information from these data. Quality control (QC) is the rst step in data analysis. FastQC is commonly used to assess data quality. We used FastQC to perform basic statistical analysis of the quality of the raw reads and the results are shown in Fig. 1.
The main purpose of data ltering is to remove low-quality data and ensure the quality of clean data. We trimmed off the sequencing adapters and low-quality fragments from the raw sequencing data using Trimmomatic software and the results are shown in Table 1.

Analysis of Clean Reads
A comparison of the clean reads with the reference genome showed that the number of unique reads accounted for ~ 50% of the total, and the ratio of repetitive sequences in the reads to total reads was ~ 10% (Table 2).  To fully detect the methylation of a genome, the NIH Roadmap Epigenomics Project (http://www.roadmapepigenomics.org/protocols) recommends a sequencing depth of 30× or more. Therefore, we used a sequencing depth of 30× to perform methylation sequencing of all samples and to calculate the coverage of each single-base site in the genome (that is, the number of reads that support that site). The coverage statistics and distribution diagrams are shown in Table 3 and Fig. 2, respectively. (3) sites_numCovg1: greater than or equal to 1× on the genome The number of bases accounted for the total length of the genome; (4) sites_numCovg5: the number of bases on the genome greater than or equal to 5 × the depth of sequencing accounted for the total length of the genome; (5) sites_numCovg10: the number of bases on the genome greater than or equal to 10 × the depth of sequencing accounted for the total length of the genome proportion.
The coverage level of the C site is an important indicator of the sequencing depth in the methylation test. We calculated the coverage of the C site separately, as well as the coverage of the C site under each context (CpG, CHH, CHG) (that is, the number of reads that support that context). The results and cumulative distribution diagram are shown in Table 4 and Fig. 3, respectively. Among the three methylation environments (CG, CHG, and CHH), CG had the highest coverage (47.99-51.59%), followed by CHG (18.73-24.75%) and CHH (2.34-4.59%).

Types and Distribution of Methylation in the high-and low-TPC lines
We performed statistical analysis of the methylated C sites in three sequence environments (CG, CHH, CHG, where H stands for A, C, or T) in the two plant parts (young buds and stem apex) for the two sets of samples: high-TPC (YC4, YC33) and low-TPC (YC11, YC15) lines (Table 5 and Fig. 4). The percentage of methylated C sites in the three sequence regions of all samples represents the percentage of the total number of C sites in that region. The percentages of CG sites were the highest (33.32-37.32%), followed by CHG (21.05-26.69%) and CHH (2.88-6.52%). The overall average level of methylation of the three sequences in all samples (from high to low) was CG (43.99%), CHG (29.09%), and CHH (26.92%).   ). An analysis of the methylation level of each chromosome in each sample showed that the ratios of methylated C sites to total methylated C sites on the chromosome in different sequence environments were consistent across all materials: CG had the highest ratio, followed by CHG and CHH (Fig. 5).

Methylation Density on Chromosomes and Distribution of Methylation Levels in Genes
B. napus is a heterotetraploid with two sets of chromosomes (A and C). Circos diagrams are commonly used to show the distribution of methylation density on chromosomes [33][34]. There were signi cant differences in the distribution of methylation density on C vs. A chromosomes in all samples. Even though A chromosomes contain more genes than C chromosomes, the methylation density of C chromosomes was generally higher than that of A chromosomes (Fig. 6).
We analyzed the average methylation levels of C sites in the CG, CHG, and CHH contexts in various functional genomic regions (such as promoter, exon, intron, 5′ UTR, 3′ UTR, and so on). Here, the 2-kb region upstream of the transcription start site (TSS) was considered to be the promoter region.
The distribution of the average methylation levels in the functional elements of the genes is shown in Fig. 7. In both the high-and low-TPC lines, the results were consistent in all contexts. The promoter regions had the highest methylation levels, followed by introns. We also analyzed the methylation levels in the upstream and downstream regions of each gene in each sample. Speci cally, we calculated the average methylation level of the C site in each gene body, 2 kb upstream of the TSS, and 2 kb downstream of the transcription termination site in each context; the results are shown in Supplementary Fig. 4. The average methylation levels were higher in the regions 2 kb upstream and downstream of the TSS than in the gene body regions in all samples.

Comparison of the Overall Methylation Levels of High-and Low-TPC Lines
We constructed a Circos diagram to display differences in the methylation levels in different samples from the high-vs. low-TPC lines [33,35]. We detected signi cant differences in the overall methylation levels of multiple chromosomes between the two groups of lines ( Supplementary Fig. 5A We also examined the differences in overall methylation levels between stem and ower buds in both groups. The C methylation level of CG showed the greatest difference between stems and young buds, followed by CHG. By contrast, for all samples, there was no signi cant difference in the level of C methylation in the CHH context between stems and young buds (Fig. 8).

DMR Analysis of High-vs. Low-TPC lines
We used DSS software to identify DMRs between the high-and low-TPC lines [36][37][38]. We generated a Circos diagram to visualize the distribution of DMRs in the genome and the results of mapping analysis. We detected numerous DMRs in all three sequence contexts (CG, CHG, and CHH) between the two sets of lines in both the stem apex and ower buds ( Supplementary Fig. 6). All comparisons generated consistent results. The most signi cant DMRs included both hypermethylated and hypomethylated regions in the CG sequence context (Supplementary Fig. 6).
We performed statistical mapping of the DMR anchoring areas (such as promoter, exon, intron, CGI, CGI shore, repeat, TSS, and TES, regions) to distinguish hypermethylated vs. hypomethylated DMRs (Fig. 9). We detected both hypermethylated and hypomethylated DMRs in all functional elements in both the high-and low-TPC lines in the stem apex and ower bud. Notably, there were many DMRs in the promoter regions of all samples, including the CG, CHG, and CHH contexts.

DMR analysis of the promoters of high-vs. low-TPC lines
Methylation of the promoters of plant genes plays an important role in regulating gene expression, thus affecting plant traits [39]. Therefore, we identi ed genes with signi cantly different methylation levels in their promoter regions. In the stem apex, 427 genes were hypermethylated and 547 genes were hypomethylated in the promoter regions in the high-TPC lines compared to the low-TPC lines ( Fig. 10C and D). In ower buds, 367 genes were hypermethylated and 531 genes were hypomethylated in the promoter regions of the high-TPC lines ( Fig. 10A and B). Compared to the low-TPC lines, 293 genes were hypermethylated and 398 genes were hypomethylated in the promoter regions in both the stem apex and ower buds (Fig. 10E). In addition, for a few genes, the promoter regions were signi cantly hypermethylated or hypomethylated (Fig. 10E).

DMRs in Promoters Underlie Differences in Gene Expression Between High-and Low-TPC Lines
Increased methylation in the promoter region of a gene inhibits gene expression, while demethylation in the promoter region of a gene promotes gene expression [39]. We combined the current results with our previously reported transcriptome data [1] and identi ed 26 genes in the high-TPC materials that were signi cantly hypermethylated in their promoter regions compared to the low-TPC materials in B. napus ower buds, resulting in signi cant decreases in their expression ( Fig. 11A and Supplementary Table 1). We also identi ed 22 genes that were hypermethylated in their promoter regions in the high-vs. low-TPC materials in the stem apex, resulting in signi cant decreases in their expression ( Fig. 11A and Supplementary Table 1). Finally, the methylation levels of the promoters of 19 and 21 genes were signi cantly reduced in the ower bud and stem apex, respectively, in high-vs. low-TPC materials, which led to signi cant increases in their expression ( Fig. 11B and Supplementary Table 1).
Therefore, the promoters of many genes are affected by hypermethylation or demethylation in the high-vs. low-TPC materials, leading to signi cant differences in gene expression. Some of these genes are related to the TPC trait. For example, in both the stem apex and ower buds of the high-TPC lines, the promoter regions of BnaC03g53050D (UBC32), BnaA05g26660D (CYSB), and BnaA10g07880D (TCP1) had low methylation levels in the high-TPC lines, and they were expressed at signi cantly higher levels in the high-TPC vs. low TPC-lines. The promoters of BnaC03g70780D (auxin associated family protein), BnaC05g36710D (myb family transcription factor), BnaAnng09670D (SMP1), BnaA09g02000D (SDH2-2), and BnaC01g12960D (NRT1.8) were hypermethylated in the high-vs. low-TPC lines, leading to signi cantly lower expression levels in both buds and the stem apex. The promoter region of BnaC09g30490D (TAF15b) was hypomethylated in the ower buds of the high-vs. low-TPC lines, leading to higher expression in this tissue.

GO and KEGG Analysis of Differentially Methylated Genes
To investigate the functions of genes whose promoter regions were differentially regulated by methylation and were differentially expressed between the high-and low-TPC lines, we performed GO and KEGG analysis of these genes (Fig. 12). The signi cantly enriched GO terms are involved in nitrogen metabolism, biosynthesis regulation, cytoskeleton composition, and so on (Fig. 12A). KEGG analysis showed that the signi cantly enriched pathways mainly involve energy-related processes such as carbon metabolism (Fig. 12B).

MiRNAs and Methylation Jointly Regulate the TPC Trait
To examine whether miRNAs and methylation are both involved in regulating the expression of genes involved in the TPC trait, were performed a joint analysis of transcriptome, methylation, and miRNA data for the high-vs. low-TPC samples. The promoter regions of 14 important genes related to growth and development were differentially methylated and their expression regulated by miRNAs, resulting in signi cant differences in gene expression between the two groups ( Table 6)  Note: In the "Sites" column,"common"indicates signi cant differential expression in shoot tips and ower buds. " ower bud" indicates signi cant differential expression in ower buds, and "stem apex" indicates signi cant differential expression in shoot tips.

Veri cation of the methylation sequencing results
The correlation of methylation levels between samples is an important indicator of the reliability of an experiment and whether the sample selection is reasonable. The closer the correlation coe cient is to 1, the higher the similarity of the methylation patterns between samples. We used the 2kb/bin sub-sequence environment to calculate the methylation level in each bin and performed Pearson correlation analysis of the data [40]. The correlations coe cients of the CG, CHG, and CHH methylation levels between the stem apex and ower buds in the four high-and low-TPC lines were very high (Fig. 13). This nding indicates that the overall differences in methylation between the two plants parts were not signi cant and that the sequencing results are reliable.
To further verify the reliability of the sequencing results, we used a DNA Bisul te Conversion Kit (Tiangen, Beijing) to process the remaining samples from the sequencing experiment (we randomly selected YC4F and YC4S). We ampli ed segments of the differentially methylated regions of 20 gene promoters by PCR, followed by sequencing after TA cloning. We compared the sequencing results with the previously generated high-throughput sequencing data. The promoter regions of all 20 genes showed the same methylation patterns in both sets of data (Supplementary Table 2).
Therefore, the results were reliable and could be used for further analysis.

Discussion
Cytosine DNA methylation is a chemical modi cation that produces 5-mC, also known as the fth base of DNA. In plants, DNA methylation occurs in three different sequence environments: symmetric CG and CHG and asymmetric CHH, where H stands for C, A, or T. These DNA methylation patterns are stably inherited through cell division. Changes in DNA methylation can occur spontaneously or be induced by genetic factors and environmental stimuli. High-throughput DNA sequencing is used for single-base resolution transcriptome and epigenome analyses, as well as genome resequencing [41]. The integration of these omics-based data has shed light on the biological roles of the epigenome [42][43]. There is considerable intra-and inter-species variation in DNA methylation patterns. The analysis of this type of data is not limited to model plant species with compact genomes, instead extending to important agronomic crops with large, complex genomes [34,[44][45][46][47]. Natural genomic variations, such as single-nucleotide mutations and structural variations, have been used at tools in plant breeding [34,47].
In the current study, we determined that the average methylation levels of all B. napus samples (from high to low) were mCG (43.99%), mCHG (29.09%), and mCHH (26.92%). A recent study indicated that a B. napus male sterile line (7365A) and a restorer line (7365B) showed completely different methylation patterns ( Fig. 14a and b) [21]. There was no signi cant difference in the methylation levels of CG, CHG, and CHH between the sterile line and the restorer line, but for both lines, the overall methylation level was mCHH > mCG > mCHG. This nding indicates that changes in methylation patterns are strongly related to the fertility of B. napus.
In the current study, all of the B. napus samples used for sequencing showed the same trend: the methylation levels of all three sequence contexts were signi cantly higher on chromosome C than on chromosome A. This result is consistent with previous studies [21,48]. The A and C chromosomes of B. napus come from Chinese cabbage (Brassica rapa) and Brassica oleracea, respectively [48][49][50]. Based on separate studies, the overall methylation level of the Brassica oleracea genome is signi cantly higher than that of Brassica rapa [51][52]. Together, these ndings further support the notion that Brassica napus evolved from the hybridization of B. rapa and B. oleracea and that the methylation patterns of parental species B. rapa and B. oleracea were stably inherited by their offspring, B. napus.
Through joint analysis of methylation and transcriptome data, we determined that the promoter regions of many important genes are hypermethylated or demethylated, which affects the differential expression of the genes between the high-and low-TPC materials. Due to hypermethylation in the promoter region of SDH2-2, the expression of this gene was reduced in both the stem apex and buds of the high-TPC lines. Indeed, increased methylation in the promoter region affects the expression of SDH2-2 in maize seeds, thereby affecting the glyoxylic acid cycle and seed germination [53]. Therefore, the differential methylation of the SDH2-2 promoter region causes signi cant differences in gene expression between the high-and low-TPC materials, thereby affecting the glyoxylic acid cycle, leading to differences in the characteristics of the two groups of materials.
NRT1.8 is hypermethylated in both the stem apex and buds of the high TPC lines, leading to signi cantly reduced gene expression. NRT1.8 affect nitrogen use e ciency in Brassica napus, and nitrogen has major effects on TPC-related traits [54][55]. The loss of methylation in the promoter regions of CYSB and UBC32 cause these genes to be expressed at high levels in the high-TPC lines. CYSB is involved in the nitrogen stress response in spinach [56]. The Arabidopsis ubiquitin-conjugating enzyme UBC32 is an ERAD element that plays an important role in brassinolide-mediated plant growth and salt stress tolerance [57]. In the ower buds of the high-TPC materials, the TAF15b promoter was hypomethylated, leading to signi cantly higher gene expression compared to the low-TPC lines. TAF15b participates in the autonomous pathway of owers, inhibits the transcription of FLOWERING LOCUS C, and regulates the development of oral organs [58][59]. Floral organ development plays an important role in silique development, which in turn is important for determining the TPC trait.
RdDM is an important epigenetic pathway that affects many traits in plants [60][61]. MicroRNA-mediated methylation is widely involved in regulating gene function in both animals and plants. In humans, miR-29 targets DNA demethylase genes, and its activity is mediated by members of the ten eleven translocation (TET) family, which may cause the dysregulation of genes involved in key cell functions, leading to lung cancer [62]. MiR-140-5p regulates T cell differentiation and reduces autoimmune encephalomyelitis by affecting CD4 + T cell metabolism and DNA methylation [63]. miR2936 and miR398-mediated DNA methylation affect respiratory energy metabolism in Arabidopsis by regulating the expression of AGO1 and AGO4 [26]. Floral organ development plays a vital role in plant reproduction and the TPC trait. The RdDM pathway mediated by microRNAs regulates the poly(A) polymerase gene PAPS1 in Arabidopsis, thus affecting both sporophyte and pollen development [28]. Therefore, miRNA-mediated DNA methylation plays an important role in the normal growth and development of animals and plants.
Here, were compared genome-wide methylation, transcriptome, and miRNA data and found that many gene promoter regions in the high-and low-TPC materials were jointly regulated by miRNAs and methylation, leading to signi cant differential gene expression. These genes are closely related to the TPC trait. For example, the BnaCnng64040D (Lipase family protein) promoter is regulated by miR159b and is hypermethylated in ower buds of the high-TPC materials, leading to signi cantly reduced expression. Lipase family protein plays an important role in regulating lipid metabolism and oil content of B. napus [64]. The promoter region of BnaC09g30490D (TAF15b), which affects organ development, is regulated by miR167c and undergoes demethylation in the high-TPC materials, leading to signi cantly higher expression in ower buds compared to the low-TPC materials. The promoter region of BnaC04g55380D (MAP65-4) is regulated by miR6029/miR9049 and is hypermethylated, which may affect the formation of the cytoskeleton in B. napus [65]. In fact, we found that together with differences in methylation, many miRNAs that are differentially expressed between high-and low-TPC lines are closely related to the formation of TPC, such as miR159, miR827, miR319, and so on [1].
Finally, GO and KEGG analysis of genes that are regulated by methylation in their promoter regions and are differentially expressed in high-vs. low-TPC materials primarily function in metabolic processes such as carbon and nitrogen metabolism. We previously reported that miRNAs are involved in the differential expression of genes between high-and low-TPC lines. GO and KEGG analysis revealed that these genes are also mainly involved in carbon and nitrogen metabolism [1]. Therefore, the results obtained in the current study are consistent with previous results. These ndings indicate that the TPC trait is simultaneously regulated by miRNAs and methylation and is closely related to carbon and nitrogen metabolism processes.
DNA methylation is an epigenetic modi cation that affects gene expression and transposable element activity. Because cytosine DNA methylation patterns are inherited through mitotic and meiotic cell division, differences in these patterns may lead to phenotypic variation. Advances in highthroughput sequencing technology have led to the generation of abundant DNA sequence data. The comprehensive analysis of genome-wide gene expression and DNA methylation patterns has revealed the underlying mechanisms and functions of DNA methylation. In addition, various associations between DNA methylation and agronomic traits have been uncovered [66]. The information obtained from such studies could be used for crop breeding based on natural epigenomic variations in the future. In addition, arti cial epigenome editing could be used as a new breeding technique to produce new crop varieties with improved agronomic traits.

Conclusions
Signi cant differences were detected in overall methylation levels between the high-and low-TPC lines in the CG, CHG, and CHH contexts in the promoters of genes in the stem apex and ower bud. In addition, 14 important genes related to growth and development were differentially regulated between the two groups due to miRNA-mediated differences in methylation levels in their promoters. These results, together with our previously generated RNA-seq and miRNA pro ling data, indicate that both methylation and miRNAs are involved in regulating the expression of genes in nitrogen-related metabolic pathways, thereby affecting the TPC trait in B. napus.  Genomic coverage map of all samples. Above: Distribution map, the abscissa is the depth of coverage, and the ordinate is the ratio of the number of base sites corresponding to the coverage to the number of bases in the whole genome; the following gure: the cumulative distribution map, the abscissa is the coverage depth, and the ordinate is greater than The ratio of the number of base sites equal to the corresponding coverage to the number of bases in the entire genome; the lines of different colors represent different samples.

Figure 3
Partial cumulative distribution map of C-site coverage. The abscissa is the coverage depth of the C site, and the ordinate is the percentage of the C site that is greater than or equal to the corresponding coverage depth as a percentage of the total C site, and the different colors represent different contexts.

Figure 4
Proportional distribution map of methylation C site. Different colors represent methylation C sites under different contexts, and the size of each part represents the proportion of methylated C sites under the corresponding contex.

Figure 6
Circos diagram of chromosome methylation density. From the outside to the inside, the CG sequence environment methylation density, CHG sequence environment methylation density, CHH sequence environment methylation density, TE original proportion density heat map, gene number density heat map; internal scale: MC density Hot label: from green to yellow to red, the methylation density is from low to high, and the TE ratio is hot: from green to black to red, the repeat sequence ratio is from low to high, and the gene density is hot: from gray to black. Indicates the number of genes from low to high. TE density calculation method: Calculate the TE (repeat original) ratio, which is the total length ratio of the repeat in the bin (the analysis is not performed without repeat); Gene density calculation: Calculate the number of genes contained in each bin.

Figure 8
Difference in methylation level between stem and ower bud. From the outside to the inside, the circle indicates the methylation level of the treatment group, the difference of methylation level between the sample groups, and the methylation level of the control group. The internal scale: the DNA methylation level indicates the level of methylation. The DNA methylation difference indicates the degree of difference in methylation levels between samples.

Figure 9
DMR anchoring area display of three sequence environment (CG, CHG, CHH). The abscissa represents the respective area categories, and the ordinate represents the number of DMRs of the hyper/hypo DMR in each area.

Figure 10
Differential gene for promoter methylation between high-and low-TPC materials. A: Genes with hypermethylation at ower bud in high TPC materials compared to low TPC materials; B: Genes with hypomethylation at ower bud in high TPC materials compared to low TPC materials; C: Genes with hypermethylation at stem apex in high TPC materials compared to low TPC materials; D: Genes with hypomethylation at stem apex in high TPC materials compared to low TPC materials. E: All high-and low-methylated genes in the promoter region at both sites of two groups of materials.

Figure 11
Methylation affects gene expression between high-and low-TPC materials in Brassica napus L. A: Hypermethylation in the promoter region results in a signi cant decrease in gene expression in the high-TPC materials compared to the low-TPC materialss; B: loss of methylation in the promoter region results in an increase in gene expression in the high-TPC materials compared to the low-TPC materials.

Figure 12
Analysis of GO and KEGG of genes regulated by promoter methylation in two sets of TPC materials Figure 13