Tissue-specific regulatory mechanism of LncRNAs and methylation in sheep adipose and muscle induced by Allium mongolicum Regel extracts

Allium mongolicum Regel (A. mongolicum) is a perennial and xerophytic Liliaceous allium plant in high altitude desert steppe and desert areas. Feeding A. mongolicum greatly reduced unpleasant mutton flavor and improves meat quality of sheep. We analyzed epigenetic regulatory mechanisms of water extracts of A. mongolicum (WEA) on sheep muscle and adipose using RNA-Seq and whole-genome Bisulfite sequencing. Feeding WEA reduced differentially expressed genes and long non-coding RNAs (lncRNAs) between two tissues but increased differentially methylation regions (DMRs). LncRNA and DMR targets were both involved in ATP binding, ubiquitin, protein kinase binding, regulation of cell proliferation, and related signaling pathways, but not unsaturated fatty acids metabolism. Besides, tissue specific targets were involved in distinct functional annotations, e.g., Golgi membrane and endoplasmic reticulum for muscle lncRNA, oxidative phosphorylation metabolism for adipose lncRNA, dsRNA binding for muscle DMRs. Epigenetic regulatory networks were also discovered to discovered essential co-regulated modules, e.g., co-regulated insulin secretion module (PDPK1, ATP1A2, CACNA1S and CAMK2D) in adipose. The results indicated that WEA induced distinct epigenetic regulation on muscle and adipose to diminish transcriptome differences between tissues, which highlights biological functions of A. mongolicum, tissue similarity and specificity, as well as regulatory mechanism of mutton odor.


Methylation characteristics distinguished by methylation context. Average Bisulfite sequencing
depth of bases was 17.70 ± 0.49, while sequencing coverages of bases (sequencing depth > 10 ×) were 66.42-78.82% (Table S3). Totally we obtained 42.56 Gb C sites and average coverage depth was 5.03 ± 0.20 (Table S4). Base numbers were 267.98 ± 12.69 (Mb), 1309.16 ± 45.64 (Mb) and 3870.56 ± 148.42 (Mb) for CG, CHG and CHH respectively (Table S4). Average methylation level of CG was 70.60% ± 1.59%, while that of CHG and CHH were only 0.27% ± 0.02% and 0.27% ± 0.02% (Table S4). Over 30.01% of C sites were methylated in CG context, which was higher than that in the whole genome context, CHG or CHH (> 1.67%, > 0.09%, > 0.16%). 92.17% ± 4.74% of methylated C sites was in CG context, while 6.52% ± 4.85% was in CHH context. For example, in normal adipose tissue, percentages of mCG, mCHG and mCHH were 94.29%, 1.26% and 4.45% (Fig. 2a). Methylation levels in CG context were higher than other context, e.g., in normal adipose tissue (Fig. 2b). Methylation sites in CHH context were mainly CAC (Fig. 2c), which was different from that in the whole genome sites (Fig. 2d). In CG context, single-sample C site methylation level was from 31.82% ± 7.80% to 82.72% ± 1.19%, and single-sample C site methylation density was from 86.12% ± 1.30% to 98.64% ± 0.29%, which were higher than that in CHH context and CHG context (Table S5). Methylation levels in CG context was consistent in eight gene functional regions among samples (Fig. 2e).  Figure S5a, S5b, S5c). These DMRs included 17,751 hyper-methylation regions and 2232 hypo-methylation regions. The number of CG hypo-methylation regions was more than the number of CG hyper-methylation regions in intron, repeat and CGI shore (Fig. 3a). However, CHH and CHG DMRs were mainly hyper-methylation, and had no difference among gene regions (Fig. 3b, Figure S6). In muscle, WEA induced 4,694 differentially expressed mRNAs in muscle ( Figure S7), including 2703 upregulated and 1991 down-regulated mRNAs related to 3092 DEGs, 1412 DElncRNAs including 687 up-regulated and 725 down-regulated lncRNAs. Similar with adipose, the overall methylation levels of muscle were higher in gene body and downstream 2 kb region than that in other gene regions ( Figure S8a). In addition, WEA decreased methylation levels in CG, CHG and CHH context (Fig. 4a-c). 5019 DMRs were induced by WEA, including 2787 hyper-methylation regions and 2232 hypo-methylation regions. Average lengths of CG DMRs and CHG DMRs were 60-70 bp (Fig. 4d, e). However, average length of CHH DMR was 110-120 bp (Fig. 4f), which was shorter than that in adipose. DMRs in muscle included more hypo-methylation regions in CG, CHG and CHH context ( Figure S8b, S8c, S8d).
WEA decreased transcriptome differences and increased methylation differences between adipose and muscle. WEA reduced transcriptional differences between muscle and adipose. In normal tissues, we discovered 107 differentially expressed mRNAs corresponding to 99 genes. Expression of 95 mRNAs were higher in muscle than in adipose ( Figure S9). After WEA treatment, DEGs number was decreased to only 3 (CACNG6; DUSP13; and 105,606,075 gene, Figure S10), indicating that WEA reduced transcriptome differences between muscle and adipose. In PCA analysis of 12 samples from the four experimental groups ( Figure S11), we discovered that in the first principal component (PC1) WEA reduced the PCA space between mutton muscle and adipose tissues, while the tissues differences were still maintained. In addition, the number of DElncRNAs was also decreased from 9 to 1. Effects of WEA on methylation between tissues were different from transcrip- www.nature.com/scientificreports/ tome. In normal tissues, average length of tissue CG DMRs was 110 bp ( Figure S12a). 9,760 DMRs were discovered, including 2339 hyper-and 7421 hypo-methylation regions. Hypo-and hyper-methylation distributions on gene region were different in CG, CHH and CHG context. The numbers of CG hypo-methylation were more than that of CG hyper-methylation in intro and repeat regions ( Figure S12b), while more CHG and CHH hypomethylation existed in all gene regions ( Figure S12c, S12d). Interestingly, compared with 2826 DMR genes between normal tissues, the number of DMR genes between treatment tissues was increased to 6226, indicating that WEA induced more DMR. WEA's effects on methylation were different to the decreasing tissue differences of transcriptome and lncRNAs. DMR genes showed distinct patterns among three methylation contexts. In normal tissues, we discovered 2,697 CG DMRs genes and 216 CHH DMRs genes, which had 90 genes in common (Fig. 5a). However, WEA increased the numbers of CG and CHH DMRs genes to 3185 and 4164 (Fig. 5b). Overall CG methylation levels in normal conditions showed higher methylation levels in adipose (Fig. 5c) and were not affected by scallions (Fig. 5d). By contrast, the CHH and CHG methylation levels of muscle were higher than that in adipose under normal condition (Fig. 5e, Fig-ure S13a), while WEA feeding induced higher CHH and CHG methylation levels in adipose (Fig. 5f, S13b).

Functional interpretation of DMRs genes.
Then, we did pathway and GO enrichment analyses for DMRs genes to investigate the underlying regulatory mechanism of increasing methylation difference between two tissues. Results indicated that the DMRs genes were mainly involved in energy metabolism, signaling pathways and cell proliferation. www.nature.com/scientificreports/ DMRs genes induced by WEA were significantly enriched in 21 GO terms, while that between normal tissues were only 7 GO terms (adjusted P < 0.01, Fig. 6a, Figure S14, Figure S15). 4 GO terms (ATP binding, intracellular signal transduction, negative regulation of cell proliferation and regulation of GTPase activity) were related to tissue difference and not affected by WEA treatment. The remaining 17 GO terms were induced by WEA treatment and were mainly about signaling transduction and signaling related responses, e.g., transmembrane receptor protein tyrosine kinase signaling pathway, Ras protein signal transduction, positive regulation of NF-kappaB transcription factor activity, transmembrane receptor protein tyrosine kinase activity, protein kinase  www.nature.com/scientificreports/ binding, phosphotyrosine residue binding, postsynaptic membrane, focal adhesion, cytoskeleton and caveola. Interestingly, ubiquitin protein ligase binding, nucleoplasm and positive regulation of cell proliferation were also identified as enriched GO annotations here, indicating their essentiality in tissue differences regulated by WEA. Methylation alterations of these GO were all hypo-methylation ( Fig. 6a, Figure S14, Figure S15), indicating WEA induced higher methylation level of these GO terms in adipose than that in muscle. In addition to tissue differences, we discovered tissue specific functional annotation of DMRs genes (see NF VS BF, NM VS BM in Fig. 6a, Figure S14, Figure S15). In adipose, enriched GO terms were also related to the signaling pathways and cell proliferation. In muscle, we found four enriched GO annotations, including intracellular signal transduction, ATP binding, response to exogenous dsRNA and double stranded RNA binding. dsRNA was only enriched in hyper DMR genes, indicating that WEA might regulate dsRNA specifically in muscle through enhancing methylation of related genes. Compared with adipose, enriched GO annotations in muscle were less than that in adipose.
Pathway enrichment indicated that the effects of WEA on hyper-, hypo-and overall methylation were about genetic information processing, environmental information processing, cellular processes and organismal systems, but not about fatty acid metabolic pathways (Fig. 6b). And the results were similar with GO, e.g., ubiquitin mediated proteolysis, Rap1 signaling pathway, MAPK signaling pathway, axon guidance, cholinergic synapse, adherens junction, focal adhesion, regulation of actin cytoskeleton, cell cycle, etc.
Functional annotation of lncRNAs' targets induced by WEA. We discovered that WEA induced 30,911 lncRNA-target pairs of 1218 lncRNAs (86.26% of total) and 2296 DEGs (74.25% of total) in muscle, while in adipose there were 24,483 pairs including 974 lncRNAs (83.18% of total) and 2035 DEGs (71.60% of total) (Methods). Results of normal tissue comparison were less, with only 359 lncRNA-target pairs including all the 9 DElncRNAs and 87 DEGs (81.31% of total).
LncRNAs targets were significantly enriched in GO terms related to energy metabolism, signaling transduction and cell proliferation, which were the same with the results of methylation, e.g., NAD binding, npBAF complex, NAD + kinase activity, pyridoxal phosphate binding, etc. Effects on signaling transduction were mainly about magnesium ion binding, protein kinase binding, growth hormone receptor binding, lipopolysaccharide receptor activity, regulation of cytokine secretion, etc. Furthermore, these targets were also significantly enriched in transcription regulation, e.g., ubiquitin protein ligase binding, ubiquitin conjugating enzyme activity, SWI/ SNF complex and transcription factor TFIID complex. In addition, lncRNAs target were also enriched in positive www.nature.com/scientificreports/ www.nature.com/scientificreports/ regulation of cell proliferation. Besides, we discovered that lncRNAs targets were correlated with macrophage activation (a biological process related to immunity). lncRNAs targets showed tissue differences in environmental information processing, cellular processes, metabolism and organismal systems pathways. For instance, the targets in muscle were significantly enriched in protein synthesis related GO terms, e.g., endoplasmic reticulum unfolded protein response, Golgi membrane, endoplasmic reticulum. These enriched GO terms were specific identified in muscle, which was distinguished from adipose. A metabolism pathway (oxidative phosphorylation pathway) was notably specific for lncRNA regulation in adipose, which was different from the results of methylation. Besides, other enriched pathways were the same with the results of methylation, e.g., ubiquitin mediated proteolysis, endocytosis, apelin signaling pathway, cGMP-PKG signaling pathway, focal adhesion, etc.
Essential targets in co-regulated module of epigenetic regulatory networks. By integrating both methylation and lncRNA results, we constructed and analyzed the regulatory networks of tissue difference or specific tissue to discover essential regulators and genes during transcriptional processes. We obtained a regulatory network of tissue difference between normal muscle and adipose (Fig. 7), including nine differentially expressed lncRNAs (Table S6) and 23 methylated target genes. Under the regulation of lncRNAs and methylation, the expressions of all the target genes were decreased in muscle. Each lncRNA regulated seven target genes on average.
Besides, we discovered that WEA induced complex networks for specific tissue on account of more regulators and more target genes. The regulatory network for muscle had 983 regulatory interactions including 517 differential lncRNAs and 104 target genes ( Figure S16). The expressions of 68 target genes were up-regulated by the regulation of both lncRNAs and methylation, while the remaining 36 genes were down-regulated. Interestingly, we uncovered that co-regulated genes in the regulatory network, that were regulated by the same lncRNAs. For instance, MAPKAP1 gene and FOXP2 gene were co-regulated genes. MAPKAP1 gene was regulated by 135 lncRNAs, while FOXP2 gene were regulated by 121 lncRNAs. 108 of MAPKAP1 gene regulators (80.00%) were the same with the regulators of FOXP2 gene (89.26%). In another example, STXBP5 gene and WDR17 gene were regulated by 46 and 40 lncRNAs respectively. 35 lncRNAs were the same, accounting for 76.09% and 87.50% individually. These co-regulated gene pairs shared by gene pairs suggested that the intake of WEA induced similar downstream regulation mechanism of lncRNAs, which could be essential regulators.
Similar with the regulatory network of muscle, the regulatory network of adipose also suggested complex regulatory mechanism ( Figure S17). This network included 346 lncRNAs and 80 target genes. 59 of the target genes were up-regulated, while the remaining 21 were down-regulated. In adipose, co-regulated modules contained more than two target genes, that was different from the co-regulated modules in muscle. In one co-regulated module, seven target genes (THSD4, SLC9A7, PDPK1, GMDS, DLG4, CDH4 and ADGRG6) were co-regulated by 18 lncRNAs, accounting for more than 31.58% of lncRNAs in regulatory network. In another example, SLC25A12, MAP2K6, REEP1, MYBPHL, CLCN1 and PHB were co-regulated by six lncRNAs. The co-regulation of target genes in these modules suggested that DEGs induced by WEA had high function correlation in adipose.
Validating key target genes in epigenetic regulatory networks by RT-PCR. We randomly selected some target genes from epigenetic regulatory network and examined their expressions by fluorescence quantitative PCR (Fig. 8) to validate high-throughput results. MYOZ1 was a key target gene regulated by eight www.nature.com/scientificreports/ lncRNAs (LNC_001888, LNC_003112, LNC_004544, LNC_007802, LNC_007807, LNC_008243, LNC_008267 and LNC_009344) and methylation. MYOZ1 expression was specific for muscle, which was reduced by WEA. PDLIM5 was regulated by methylation and four lncRNAs (LNC_001888, LNC_004544, LNC_008267 and LNC_009344). PDLIM5 was expressed both in muscle and adipose. However, WEA effects on PDLIM5 expression were distinct in the two tissues. WEA increased PDLIM5 expression in adipose, while decreased PDLIM5 expression in muscle, suggesting different biological functions of PDLIM5 in the two tissues. Results of LMO7 were similar with PDLIM5, revealing potential functional correlations between LMO7 and PDLIM5. Expression patterns of TFRC were different from the above three target genes. TFRC expressions were reduced by WEA, which was the same in two tissues, indicating that biological functions of TFRC were not specific to muscle or adipose. Expressions of the four target genes in fluorescence quantitative PCR were in accordance with RNA-seq results, suggesting possible distinct biological functions of target genes in epigenetic regulatory networks.

Discussion
In this study, we analyzed epigenetic effects of WEA on muscle and adipose tissues with both methylation and transcriptome data to uncover how this plant reduce the unpleasant flavor of mutton. The results provided an overview of the whole genome changes induced by WEA and highlighted the underlying regulatory mechanisms of methylation and lncRNAs.
Here, we discovered that WEA reduced DEGs number of normal tissues to only 3 (CACNG6, DUSP13 and 105,606,075), suggesting that transcription differences of genes were minimized to maintain tissue specificity. These three DEGs induced by WEA might be essential genes in remaining tissue specificity. CACNG1 and DUSP13 were reported to be closely correlated with muscle in animals. CACNG1 (Voltage-dependent calcium channel gamma-1 subunit) encodes a muscle-specific isoform of the Ca 2+ channel gamma subunit. CACNG1 is increased in the lamina propria of the old rats, and the deleterious mutation is discovered in a young patient with metastatic urothelial carcinoma 33,34 . DUSP13 (Dual specificity protein phosphatase 13) encodes two atypical DUSPs, DUSP13B/TMDP and DUSP13A/MDSP, that are specifically expressed in testis and skeletal muscle respectively 35,36 . In addition to the known biological function of CACNG1 and DUSP13, our results provide more evidence for further understanding genes functions to maintain tissue characteristics under the regulation of WEA in sheep.
We discovered that the regulatory changes of lncRNAs and methylation had different patterns: decreased differences in lncRNAs and increasing differences in methylation, indicating that WEA had distinct influences on lncRNAs and methylation regulation. Targets of lncRNAs and methylation were both significantly enriched in energy metabolism, signaling pathways and cell proliferation related annotations, e.g., ATP binding, ubiquitin protein ligase binding, protein kinase binding, positive regulation of cell proliferation, focal adhesion, etc. These annotations were closely correlated with muscle and adipose physiology. For instance, focal adhesion is a dynamic cellular process in cell proliferation and migration, and the phosphorylation of focal adhesion molecules is significantly reduced by talin modulator to inhibit proliferation of human aortic vascular smooth muscle cells 37 .
Besides, these two types of epigenetic regulations also had differences. On one hand, we discovered that either hyper-, hypo-or overall methylation were not related to metabolic pathways, indicating that methylation induced by WEA did not regulate metabolism directly. In addition, it was notable to uncover that WEA induced higher methylation of dsRNA related genes in muscle than that in adipose, indicating functional correlations between dsRNAs and methylation. dsRNAs can be cleaved into small interfering RNAs (siRNA) to trigger RNA interference in plants, animals and some fungi, and is related to methylation. For instance, the formation of dsRNAs for Adar1 binding to Socs3 mRNA is influenced by Tet2 deficiency probably through cytosine methylation-specific readers in mammals 38 . Also, H3 mK9 is the first epigenetic modification made in response to dsRNA-derived species, with DNA methylation accumulating later in the process 39 . Therefore, these dsRNA related DEGs (TLR3, MAPK1, DHX15, VIM and TFRC) induced by WEA might be essential targets that could also induce downstream DNA methylation and be further investigated to uncover the molecular mechanism in www.nature.com/scientificreports/ detail. Other methylation specific targets were also enriched in the regulation of GTPase activity, Ras protein signal transduction, positive regulation of NF-kappaB transcription factor activity, transmembrane receptor protein tyrosine kinase activity, etc. On the other hand, the targets of lncRNAs were specifically enriched in magnesium ion binding, growth hormone receptor binding, lipopolysaccharide receptor activity, transcription factor TFIID complex, etc. Besides, the function of targets of lncRNAs also suggested tissue specificity. Targets were significantly enriched in protein synthesis related GO terms for muscle, e.g., Golgi membrane, endoplasmic reticulum. One metabolism pathway of oxidative phosphorylation was notably specific for adipose. These annotations were distinct between the regulation of lncRNAs in the two tissues.
Another highlight of this work was combining lncRNAs and methylation regulation together to construct the essential epigenetic regulatory network of WEA and tissue differences. In the epigenetic regulatory network of tissue difference, the expressions of all the target genes were lower in muscle than that in adipose, e.g., ATP1A2, CACNA1S, CAMK2D, ATP2A1 and RYR1. ATP1A2 (Na, K-ATPase α2 isoform) is abundantly expressed in skeletal muscle and is regulated by muscle use to maintain contraction and resist fatigue of muscles 40 . ATP1A2 is mapped to meat quality quantitative trait loci (rs344748241) in a Duroc pig population, demonstrating that ATP1A2 is highly associated with muscle electric conductivity 41 . Thus, these genes were identified as essential genes to maintain tissue differences under normal conditions by lower expression in these signaling pathways.
As for the epigenetic regulatory network of each tissue, the targets in co-regulated modules all played an important role in signaling transduction, transcription regulation and tissue development. For example, in the co-regulated module of muscle, MAPKAP1 was such an essential gene that were co-regulated by both lncRNAs and methylation to affect its downstream biological processes, e.g., mTOR signaling pathway, Ras protein signal transduction and regulation of transcription. FOXP2 was co-regulated with MAPKAP1, and was involved in transcription, skeletal muscle tissue development and smooth muscle tissue development. In addition to the importance of these co-regulated targets, our results indicated that these two key genes were regulated by the lncRNAs and methylation regulation that were induced by WEA. Also, MAPKAP1 and FOXP2 might had closely functional correlations to reduce meat quality for sheep. Except common signaling pathways, we also noticed that these co-regulated genes (PDPK1, ATP1A2, CACNA1S and CAMK2D) were related with Insulin secretion, resistance and signaling pathways, indicating that WEA would also influence insulin to activate downstream processes. Co-regulation of target genes in these modules suggested that WEA induced highly functional correlation of DEGs.

Conclusions
This study provides a comprehensive overview of biological responses induced by WEA in sheep muscle and adipose. Results suggested that intake of WEA did not directly regulated lipid metabolism but regulated changes of lncRNAs and methylation to affect energy metabolism, signaling pathways and cell proliferation related processes, which were upstream regulating pathways of lipid metabolism. Therefore, this work highlights the underlying molecular mechanism of biological function of A. mongolicum, tissue specificity and similarity, as well as the regulation of mutton odor and animal nutrition.

Preparation of water extracts from A. mongolicum. A. mongolicum powder (Haohai Biotechnology,
Inner Mongolia, China) was de-esterified by petroleum ether for 6 h at room temperature. After drying at 65 °C, it was mixed with distilled water at a ratio of 20 ml/g and oscillated by constant temperature water bath shake at 80 °C for 8 h. After filtration, residue was discarded, and solutions were dried in constant temperature drying oven at 65 °C to concentrate to 1/3 of the original volume. Samples were added by anhydrous ethanol of the same volume, kept at 4 °C overnight, and centrifuged at 4000 rpm for 30 min to collect the precipitation. After freezedrying, water extracts of A. mongolicum (WEA) were obtained.
Feeding and management of experimental animals. All the management and experiments of animals used in the study were approved by Institutional Animal Care and Use Committee of Inner Mongolia Agricultural University and were in accordance with the ARRIVE Guidelines 42 . Experiments were conducted at Fuchuan Feed Science and Technology Co., Ltd. (Bayannaoer, Inner Mongolia, China). Thirty healthy, female, Duhan hybrid meat sheep (4.5-month-old, 36 ± 3.5 kg body weight) were randomly divided into control (B) group and treatment (N) group equally. Sheep were fed with basal diets and mixture of basal diets and WEA (3.4 g per day per animal). Experiments lasted for 75 days (from September 10 to November 24, 2017) including a 15-day preliminary feeding period for adaptation and a 60-day experimental feeding period. Experimental sheep were fed twice a day. Water and food were available ad libitum. Compositions of basal diet were listed in Table S7, and nutrient levels of basal diet were listed in Table S8.
After 24 h pre-slaughter fasting time, 6 animals were randomly selected from B and N groups and sent to the slaughterhouse. Longissimus dorsi muscle (M) and subcutaneous fat (F) were immediately taken after slaughter and washed with normal saline. Tissue samples (BM, BF, NM and NF) were cut into small pieces and packed into multiple cryopreservation tubes respectively and stored in liquid nitrogen and transported by dry ice.
RNA isolation, library preparation and sequencing. Total 50 . Transcripts predicted with coding potential by the four software were filtered out, and the rest without coding potential were candidate lncRNAs.

Quantification of expressions and computation of differential expression. Cuffdiff in Cufflinks 45
was used to calculate FPKMs of lncRNAs and coding genes. Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group. Ballgown suite includes functions for interactive exploration of the transcriptome assembly, visualization of transcript structures and feature-specific abundances for each locus, and post-hoc annotation of assembled features to annotated features 51 . Transcripts with P-adjust < 0.05 were defined as differentially expressed. Expression profiles of DEGs were illustrated by Heatmap using pheatmap R package 52 .
Prediction of target genes of lncRNAs. Potential target genes of lncRNAs were predicted by co-location and co-expression. For co-location, target genes were defined as the genes located within the region between 100 kb up-and down-stream of lncRNAs. For co-expression, lncRNAs were co-expressed with their target genes. Expression correlation of lncRNAs and genes was computed and Pearson correlation co-efficient (PCC) was applied. P-values were adjusted by Benjamini and Hochberg false discovery rate. Co-expressed lncRNAs and genes were defined as the pairs with |PCC|> 0.95 and adjust P value < 0.05. Potential target genes were defined as genes satisfied with either one of two criterions.
Library preparation and quantification for Bisulfite sequencing. Total DNA was extracted from muscle and adipose using Trizol (TaKaRa Biotechnology Co., Ltd.). Genomic DNA purity and concentration was measured using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and Qubit DNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA) respectively. Bisulfite sequencing was performed as published procedures 53 , 5.2 µg genomic DNA spiked with 26 ng lambda DNA were fragmented by sonication to 200-300 bp with Covaris S220, followed by end repair and adenylation. Cytosine-methylated barcodes were ligated to sonicated DNA. DNA fragments were treated twice with bisulfite using EZ DNA Methylation-GoldTM Kit (Zymo Research), before resulting single-strand DNA fragments were PCR amplificated using KAPA HiFi HotStart Uracil and ReadyMix (2X). Library concentration was quantified, and insert size was assayed on Agilent Bioanalyzer 2100 system. Library preparations were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated.
Processing of methylation sequencing data. FastQC (v0.11.5) 54 was used to perform basic statistics on read quality. Reads sequences in FASTQ format were pre-processed through Trimmomatic (v0.36) 55 . Remaining reads were counted as clean reads for all subsequent analyses. Bismark software (v0.16.3) 56 was used to perform alignments to a reference genome. Reference genome was transformed into bisulfite-converted version and indexed using bowtie2 57 . Reads were transformed into fully bisulfite-converted versions. Sequence reads that produced a unique best alignment from the two alignment processes were compared to genomic sequence and methylation state of all cytosine positions in the reads was inferred. Sequencing depth and coverage were summarized using deduplicated reads. Sodium bisulfite non-coversion rate was calculated as the percentage of cytosine sequenced at cytosine reference positions.

Estimation of methylation levels and differentially methylated regions (DMRs).
To identify methylation site, we modeled sum Mc of methylated counts as a binomial (Bin) random variable with methylation rate r: We divided sequence into 10 kb bins. and calculated sum of methylated and unmethylated reads counts in each window. Methylation level (ML) for each window or C site was fraction of methylated Cs: mC ∼ Bln(mC + umC * r ML(C) = reads(mC) reads(mC) + reads(C)