Integrated Study of Transcriptome-wide m6A Methylome Conrms a Potential Mechanism for Transcriptional Regulation During Yak Adipocytes Differentiation

Background: Yak (Bos grunniens) is considered as an iconic symbol of Tibet and high altitude, but they suffer from malnutrition during the cold season that challenges the metabolism of energy. Adipocytes perform crucial role in maintaining the energy balance and adipocyte differentiation is a complex process involving multiple changes in expression of genes. N 6 -methyladenosine (m 6 A) plays an dynamic role in post-transcription gene expression regulation as the most widespread mRNA modication of the higher eukaryotes. However, currently there is no research existing on the m 6 A transcriptome-wide map of bovine animals and their potential biological functions in adipocyte differentiation. Results: We performed methylated RNA immunoprecipitation sequence (MeRIP-seq) and RNA sequence (RNA-seq) to determine the distinctions in m 6 A methylation and gene expression during yak adipocyte differentiation. In yak adipocyte and preadipocyte the content of m 6 A and m 6 A associated enzymes were substantially different. In the two groups, a total of 14,710 m 6 A peaks and 13,388 m 6 A peaks were identied. In the most part, m 6 A peaks were enriched in stop codons, 3’-untranslated regions and coding regions. The functional enrichment exploration displayed that differentially methylated genes participated in some of the pathways associated with adipogenic metabolism. In addition to that there was a positive association between m 6 A abundance and levels of gene expression, which displayed that m 6 A plays vital role in modulating gene expression during yak adipocyte differentiation. Conclusions: In short, it could be concluded that the current study provides a comprehensive explanation of the m 6 A features in the yak transcriptome, offering in-depth insights into m 6 A topology and associated molecular mechanisms underlying bovine adipocyte differentiation which might be helpful for further understanding its mechanisms.

of the m 6 A features in the yak transcriptome, offering in-depth insights into m 6 A topology and associated molecular mechanisms underlying bovine adipocyte differentiation which might be helpful for further understanding its mechanisms.

Background
The yak is the major livestock breed on the Qinghai-Tibet Plateau, and is the only large ruminant domestic specie that can survive in harsh environments along with high altitude (3,000 ~ 5,500 m), severe cold, strong UV radiation and in low oxygen [1,2]. Because of outstanding adaptability and production output, yak performs a key role in the Tibetan region. Approximately, 15 million or more than 90% of the world's overall yak population are living in Chinese territory, ensuring daily necessities for local herders such as meat, milk, wool, skins, fuel and economic bene ts [3,4]. In the Qinghai-Tibet plateau, domestic yaks mainly grow on the natural herbage of native pastures below typical grazing conditions [5]. Owing to the seasonal variability of herbage, the yak must constantly undergo insu cient feeding during the harsh winter season (October-May), which leading to the large variations in seasonal body weight and circular rhythm 'live in summer, weighty in autumn, thin in winter and dead in spring' [6]. Consequently, the subcutaneous adipose layer of yak accumulates rapidly in the summer and early autumn, providing essential energy requirements and withstanding severe cold by selective fat catabolism in cold season [7]. The adaptive mechanism allows the yak as a fascinating model for research into the adipose metabolism of plateau domestic animals. N 6 -methyladenosine (m 6 A) was rst discovered in the 1970s as the most prevalent internal modi cation of polyadenylated mRNAs and long non-coding RNAs (lncRNAs) in higher eukaryotes [8][9][10][11][12][13]. The modi cation of m 6 A methylation is mounted by a series of m 6 A methyltransferases labeled as writers: methyltransferase like 3 and 14 (METTL3 and METTL14), Wilms Tumor 1-associated protein (WTAP), VIRMA, vir-Like m 6 A methyltransferase associated (KIAA1429), RNA Binding Motif Protein 15 (RBM15), and zinc nger CCCH domain 13 (ZC3H13) [14][15][16][17][18][19][20][21]. Beside, m 6 A demethylases eliminate methylation from RNAs to enable a delicately dynamic equilibrium modi cation, and named as eraser: fat mass and obesity-associated protein (FTO) and α-ketoglutarate-dependent dioxygenase alkB homolog 5 (ALKBH5) [22,23]. Besides, a category of proteins called readers that can identify information of RNA methylation modi cation and engage in downstream mRNA translation, degradation, binding microRNA and RNAprotein interaction [24][25][26][27]. The speci c proteins primarily contain the YTH domain family (YTHDF1- 3) and IGF2BPs (IGF2BP1-3) [28][29][30][31]. Currently, m 6 A modi cations have been reported in several area of RNA metabolism, such as RNA localization, transport, splicing, stability and translation [24][25][26][27].
Moreover, previous studies described that m 6 A modi cation of mRNA performs a signi cant biological function in controlling cell metabolic processes. It dominates the transition fate of cells in mammalian embryonic stem [32], regulates the differentiation and initiation of meiosis in murine spermatogonial stem cells [33], and maintains the myogenic potential in proliferative skeletal muscle progenitors [34]. In particular, FTO facilitates the differentiation of mouse pre-adipocytes by regulating alternative splicing of pre-mRNAs with adipogenesis-related genes [35]. As well as, a recent study revealed that RNA m 6 A modi cation has a potential function in the deposition of porcine adipose tissue [36]. Considering the notable functions of m 6 A modi cation mentioned above, it is logical to assume that m 6 A modi cation may also refer to bovine adipocyte differentiation. Previously, there is no any documented relationship lies between m 6 A modi cation and bovine adipocyte differentiation. Therefore, it is necessary to detect m 6 A sites on transcriptome-wide level for identi cation the potential biological function of RNA m 6 A modi cation. In 2012, two independent researches established a method of m 6 A RNA immunoprecipitation accompanied with high-throughput sequencing (MeRIP-seq), followed the rst discovery of the N 6 -methyladenosine modi cation map to methylomes of m 6 A RNA with a resolution of 100-nucleotides [37,38]. Simultaneously, MeRIP-seq has been used to identify the m 6 A pro le in humans and mice. These results revealed that m 6 A is predominantly located close to stop codons, 3'-untranslated regions (3'-UTRs), and also in long internal exons and transcription start sites, indicating it performs a crucial role in the regulation of gene expression after transcription. These innovative researches proved that the construction of transcriptome-wide m 6 A methylome maps is suitable for further exploring the possible biological modi cation features in a new area of research.
In order to further explore the role of m 6 A modi cation in bovine adipose metabolism, and to facilitate m 6 A studies in plateau domestic livestock, we initially extracted preadipocytes from yak adipose tissue and induced them into mature adipocyte. Subsequently, the rst known transcriptome-wide m 6 A methylome pro le was obtained in yak with MeRIP-seq, and the m 6 A modi cation features were elucidated during yak adipocyte differentiation. Differential m 6 A RNA modi cations between yak preadipocyte and mature adipocyte were veri ed that regulate gene expressions and pathways related to energy metabolism of adipose. For our research, it can be conjectured that the m 6 A modi cation could be a breakthrough point for exploring the yak energy metabolism.

Results
The yak preadipocyte induced differentiation and global m 6 A quanti cation The result of Oil Red O showed that the visibility of lipid droplets in adipocytes increased signi cantly at day 12 compared to day 0 after induction with adipogenic agents (Fig. 1a-b). Meanwhile, the expression of adipocyte differentiation speci c marker genes (PPARγ, C/EBPα and FABP4) was signi cantly elevated on day 12 (adipocyte) than day 0 (preadipocyte) (Fig. 1c), suggesting preadipcyte full differentiation into adipocyte. Subsequently, to overview the m 6 A methylation during yak adipocyte differentiation, the expression of RNA methylation-related genes were contrasted by quantitative real-time PCR (qRT-PCR) detected, including METTL3, WTAP, METTL14, FTO, ALKBH5 and YTHDF1/2/3. Comparing the group of preadipocytes (Pread0) and adipocytes (Ad), the ndings showed that the expression level of methyltransferases (METTL14, WTAP and METTL3) and ALKBH5 were dramatically upregulated, while FTO was substantially downregulated and m 6 A-binding proteins (YTHDF2 and YTHDF3) was drastically upregulated (Fig. 1d). Furthermore, the content of m 6 A in the group of adipocyte was signi cantly higher compared with the preadipocyte group (Fig. 1e). Thereby, we hypothesized that during yak adipocyte differentiation the difference of m 6 A methylation may exist, in which furtherly discovered using MeRIPseq.
Transcriptome-wide m 6 A-Seq revealed global m 6 A modi cation patterns during yak adipocyte differentiation The yak adipocyte and preadipocyte of three biological replicates were used for transcriptome wide m 6 Asequencing (m 6 A-seq) and RNA-sequencing (RNA-seq) assays. In total, 12 libraries were sequenced, comprising of three replicates of preadipocyte and adipocyte for input and MeRIP samples (Additional le 1). With each MeRIP library, an average of 9.22 giga base-pair (Gb) of high-quality data was produced, and 9.49 Gb per input library (RNA-seq dataset). Then, we eliminated reads containing adapter pollutants, low-quality and indeterminate bases, an average of 7.17 Gb and 7.11 Gb obtained from per MeRIP and input libraries, respectively. The valid data were mapped to Bos_mutus genome (Version: BosGru_v2.0) using HISAT2. The proportions of mapped reads ranged from 87.96 to 96.57%, correspondingly (Additional le 1). In the yak Ad group, R package exomePeak found a total of 14,710 m 6 A peaks, containing transcripts of 9633 genes. Likewise, 13,388 m 6 A peaks were found in the Pread0 group corresponding to transcripts of 9142 genes ( Fig. 2a-b). In addition, 5848 peaks were consistently observed in two groups, and 3964 genes within the groups were modi ed by m 6 A. Compared to the Pread0 group, the Ad group had 9226 new peaks occurring with the absence of 7904 peaks, re ecting the signi cant difference between Pread0 and Ad groups in global m 6 A modi cation trends (Fig. 2a-c). m 6 A methylomes were ulteriorly mapped by HOMER software to de ne whether RRACH motif (R represents purine, A is m 6 A and H is U, A or C) were ubiquitous in our detected m 6 A. The results of the enrichment analysis in both groups showed that m 6 A RRACH consensus sequences (Fig. 2d) accorded with previous studies, in which strengthens the credibility of the m 6 A peaks and con rms the presence of a prevailing methylated modi cation mechanism.

Analysis of m 6 A modi cation distribution in yak transcriptome
We analyzed metagene models of m 6 A peaks in the global transcriptome to identify the differential distribution of m 6 A in transcripts. Our ndings indicated that m 6 A peaks were strongly associated with two different co-ordinates: enriched instantly the end of the 5' untranslated regions (5' UTRs), coding sequence (CDS), stop codon and starting of the 3' untranslated region (3' UTRs) in Ad and Pread0 (Fig.   3a); whereas, the end peaks were noticeable abundant than beginning in CDS. Subsequently, to systematically calculate the enrichment, we investigated no overlapping transcript segments per m 6 A peak with 5' UTR, CDS and 3' UTR ( Fig. 3b), in which most of them were abundant in CDS. Afterwards, we explored the distribution of m 6 A modi ed peaks with each gene, nding that almost of 60% affected genes hold only one m 6 A peak and most genes conatain one to three m 6 A peaks (Fig. 3c).

Analysis of the GO and KEGG pathways of differentially methylated genes
Comparison was performed for the abundance of m 6 A peaks between preadipocytes and adipocytes. These ndings exposed that 118 markedly hypermethylated m 6 A peaks and 51 substantially hypomethylated peaks were obtained (|log2 (fold change)| > 1, p< 0.05) (Fig. 4a). The residual peaks of the m 6 A were viewed as unaltered peaks. According to Integrative Genomics Viewer (IGV) software, the differentially methylated sites in these two groups showed altered intensity, with the GGACU motif associated m 6 A peaks (Fig. 4b). Moreover, differentially methylated m 6 A peaks represented genes were investigated by GO and KEGG pathway analysis, revealing the biological signi cance of m 6 A methylation during yak adipocyte differentiation. GO analysis revealed that differentially methylated genes were mainly implicated with DNA-templated and regulation of transcription by RNA polymerase II (ontology: biological process), cytoplasm, nucleus and integral component of membrane (ontology: cellular component), and transcription factor and microtubule binding (ontology: molecular function) (Fig. 4c, Additional le 2). Meanwhile, the biological enrichment of KEGG indicated that the genes differently methylated were substantially related to the adipogenic metabolism regulation pathways, NOD-like receptor signaling pathway, FoxO signaling pathway, Ether lipid metabolism, cAMP signaling pathway and Hippo signaling pathway (Fig. 4d, Additional le 3).

RNA-seq identi cation of genes differentially expressed in both groups
An analysis of the RNA-seq dataset (m 6 A-seq input library) displayed that the trends of global mRNA expression between preadipocyte and adipocyte were considerably different. There were 648 signi cantly different mRNAs, including 300 up-regulated and 348 down-regulated (|log2 (fold change)| > 1, p<0.05), as shown in Fig. 5a. Then, we conducted a clustered heat map to further explore the potential roles of the genes (Fig. 5b, Additional le 4). Furthermore, GO ontology and KEGG pathway were performed to analyze the differentially expressed genes. As Fig. 5c and Additional le 5 displayed, the top 20 most notable functional annotations included regulation of glucose metabolic process, canonical Wnt signaling pathway, positive regulation of cell proliferation and insulin-like growth factor ternary complex, which in uenced adipocyte differentiation. Meanwhile, the pathway exploration revealed that signaling pathways regulating pluripotency of stem cells, ECM-receptor interaction, PI3K-Akt signaling pathway and FoxO signaling pathway were signi cantly enriched (Fig. 5d, Additional le 6), revealing that differentially expressed genes were potentially participated in adipogenic metabolism.

Conjoint analysis of RIP-seq and RNA-seq data with both groups
We found an interesting relationship of differentially methylated m 6 A peaks and gene expression patterns in preadipocytes and adipocytes through cross-analysis the m 6 A-seq and RNA-seq results, in which a positive correlation existed in differentially methylated m 6 A peaks and gene expression level (Fig. 6a). Otherwise, all genes were segregated mainly four types: 8 hypermethylated and upregulated genes termed 'hyper-up'; 7 hypomethylated and downregulated genes termed 'hypo-down'; 12 hypermethylated and whilst downregulated genes termed 'hyper-down'; and 2 hypomethylated and whilst upregulated genes termed 'hypo-up' (Fig. 6b), There were slightly more 'hyper-up' and 'hypo-down' than 'hyper-down' and 'hypo-up'. Table 1 lists the expression of genes were signi cantly differently (|log2 (fold change)| > 1, p<0.05), comprising signi cantly differently methylated peaks. Then, both groups were evaluated the overall expression levels of the m 6 A-methylated and non-m 6 A-methylated transcripts (Fig. 6c). These suggested that in yak adipocyte differentiation, m 6 A modi cations appear to have a positive association of mRNA expression. Further, we were wondering if the position of m 6 A peaks on RNA transcripts or the number of m 6 A peaks per transcript is correlated with the levels of gene expression. Based on m 6 A modi cation sites, RNA transcripts were classi ed into subgroups. As shown in Fig. 6d, m 6 A modi cations of RNA transcripts in CDS, 5'UTR or 3'UTR do not clearly difference with gene expression.
Through studying m 6 A-modi ed sites and relative expression levels of genes, revealing that the genes have three or four modi ed sites appears to be more abundant in contrast with other m 6 A-modi ed sites (Fig. 6e). Table 1 List of 29 genes with signi cant changes in m 6 A and mRNA transcript abundance in yak adipocyte as compared to preadipocyte. We implemented qRT-PCR to con rm the expression of differentially methylated genes between adipocyte and preadipocyte. The mRNA expression pattern was consistent with the RNA-seq data ( Fig. 7a-b), which con rms the validity of our transcriptome results.

Discussion
The harsh environment of Qinghai-Tibet plateau encourages yak to develop a special mechanism for energy metabolism. As an organ for energy metabolism, adipose tissue plays an crucial role in this process. To date, it has been found that epigenetic regulation is engaged in various biological processes including embryo development, stem cell self-renewal, DNA damage response, primary miRNA processing and energy metabolism [39][40][41][42][43]. In recent years, as the most extensive and plentiful internal modi cation on mRNAs, m 6 A modi cation has been a major focus in the area of epigenetic regulation [44].
Furthermore, the potential roles of m 6 A modi cation in most domestic animals, and especially for adipogenic differentiation, remained largely unknown. For the rst time, our study established a comprehensive transcriptome-wide pattern of m 6 A modi cation in yak preadipocyte and adipocyte using MeRIP-Seq technology to explore the function of m 6 A modi cation in bovine adipogenic differentiation.
Our ndings showed that yak mRNA m 6 A sites were primarily located in CDS, 5'UTRs and 3'UTRs and the distribution semblable with humans and mice [28,45], suggesting that in mammalian transcriptomes the overall distribution of m 6 A sites is in similar manner. According to previous studies, the consitent motif pattern of "RRACH" was over-represented the m 6 A motif sequence area [38,46,47]. Accordingly, in comparison with previous studies [37,38], the consensus motif (GGACU) sequence in the yak transcriptome was appropriately identi ed revealing that RNA adenosine methylation was conserved in mammals.
Earlier studies indicated that m 6 A modi cation is closely related to gene expression [45,[48][49][50]. GO analysis explored the differentially methylated genes, which participated in the transcript regulation with a variety of transcription factors by RNA polymerase II. For example, FOXO1 identi ed as a Forkhead transcription factor controlling the differentiation of adipocytes [52], and many members of the ZNF family considered as the crucial eukaryotic transcription factors involved in adipogenic metabolism [53], indicating m 6 A methylation participate in lipid metabolism. The KEGG pathway analysis revealed that the signaling pathway of differentially methylated genes closely related to adipose metabolism, such as FoxO signaling pathway, Ether lipid metabolism, Glycerophospholipid metabolism and Hippo signaling pathway-multiple species. Especially, FOXO1 furtherly found involved in FoxO signaling pathway which demonstrated important for adipocyte differentiation [52]. As a TEA domain family transcription factor, TEAD4 selected from Hippo signaling pathway-multiple species, in which recruits the cofactors VGLL4 and CtBP2 to inhibit murine adipogenesis [54]. To summarize the above ndings, we concluded that activating the FoxO signaling pathway and Hippo signaling pathway through m 6 A methylated gene may perform a key function during differentiation of yak adipocyte.
Integrated analysis of m 6 A-seq and mRNA-seq data exposed that 29 signi cant change genes existed in adipocyte group with differently methylated m 6 A sites compared to preadipocyte. Several of the genes have been con rmed to regulate adipose metabolism and adipogenic differentiation, such as ZNF395, KLF9, TEAD4, FOXO1, and UHRF1. ZNF395, the mRNA of which is hypermethylated and the expression upregulated in adipocyte group compared to the preadipocyte group. As a member of the C2H-type Zinc nger proteins, ZNF395 is classi ed to Papillomavirus-binding factor and Huntington disease gene regulatory regionbinding protein 2 [55]. Experiments of loss and gain function demonstrated that ZNF395 interacts with PPARG2 to modulate the transcriptional regulatory pathway that may be necessary for preadipocyte differentiation [56]. Besides, previous literature reported that mesenchymal stem cells were co-transduced with ZNF395 and PPARG2 enhanced the endogenous expression of PPARG2 and C/EBPα which are necessary for adipocyte differentiation [56,57]. In addition to that it was reported that Krüppellike factor 9 (KLF9), deemed to the basic transcription element-binding protein-1 (BTEB1) could transactivate PPARγ2 to regulate adipogenesis in the 3T3-L1 cell line [58]. Besides, Kimura Hiroko et al. found that KLF9 triggered the early stage of adipogenesis by promoting the C/EBPβ gene expression in 3T3-L1 cells [59]. Ubiquitin-like with PHD and RING Finger domains 1 (UHRF1) has been widely documented to promote cell proliferation [60]. Additionally, a study revealed that UHRF1 facilitates the proliferation of human adipose-derived stem cells and represses adipogenesis via inhibiting peroxisome proliferator-activated receptor γ [61]. These ndings suggested that m 6 A modi cations may perform an essential role during yak adipocyte differentiation.

Conclusions
Current ndings displyed that the m 6 A distribution patterns in the yak transcriptome and close association between m 6 A methylation and adipose metabolism. Furthermore, we have also explore the correlation between the m 6 A methylation and the level of gene expression, indicating a regulatory role for m 6 A in adipocyte differentiation. These results provides additional knowledge of m 6 A methylation in adipose tissues and it set the foundation for further understanding its possible roles and regulatory mechanisms, which could be helpful for exploration the yak adaptive mechanism in harsh environment.

Pre-Adipocyte Isolation
The Datong Yak Breeding Center (Datong County, Qinghai, China) provided three healthy 3-day-old Datong yaks. The night before slaughter, yaks were not fed. On the next morning, the yaks were humanely sacri ced by the way of electrical stunned (90 V, 10 s, and 50 Hz) at a commercial slaughter facility and exsanguinated for necessary to ameliorate the suffering, according standard approved industry protocols. The subcutaneous adipose tissue was harvested according to the protocols and guidelines for animal ethics of the People's Republic of China. Yak adipocytes were isolated from subcutaneous fat tissue that was ushed with penicillin (200 U/mL) and streptomycin (200 U/mL) added to the phosphate saline buffer (HyClone, Thermo Fisher Scienti c, Carlsbad, CA, USA). After that they were nely minced into about 1 mm 3 pieces with an aseptic setting. The segments were digested by Type I collagenase in a continuously agitated water bath at 37°C for 60-90 min. With a 40μm nylon mesh lm, indigestible material was screened and the lter liquor was resuspended for 5 min at 1400g. The sediment was subsequently hatched at room temperature for 10 min with the erythrocyte lysis buffer (0.154 M NH4Cl, 10 mM KHCO 3 , 0.1 mM EDTA). The cells were then ltered with 200μm nylon mesh lm, and rinsed twice with a serum-free medium. After 5 minutes of centrifugation at 1400g, preadipocytes were harvested and solubilized in the growth media, including DMEM-F12 (Hyclone, USA) supplemented with 10% fetal bovine serum (FBS, Gibco, USA).

Adipogenic differentiation and staining of Oil red O
To induce adipogenic differentiation, pre-adipocyte was induced for 2 days by adipogenic compounds composed with 3-isobuty-methylxanthine (MIX) (Sigma, USA), dexamethasone (Sigma, USA), rosiglitazone (Sigma, USA), and insulin (Sigma, USA) after cells con uence approached 70% in growth media. The medium was replaced after two days with DMEM-F12 containing 10% FBS, penicillin (200 U/mL), streptomycin (200 U/mL) and 5 ng/mL of insulin, and updated with cycles of two days, until day 12. The cells were usually ushed twice with PBS and set for 1 h in 4% formalin. Cells were then reacted at room temperature for 30 min with a saturated solution Oil Red O. And then, cells were rinsed three times with sterile water, and photographs were acquired from light microscopy.

Quantitative real-time PCR
Real-time RT-PCR was accomplished in a CFX Link Real-Time PCR Detection System, and 10 μL volume of reaction consisting of 5 μL 2xSYBR Premix Ex Taq II, 0.4 μL primers (10 μM), and 0.8 μL cDNA. The reaction condition was as follows: denaturation for 30s at 95°C followed by 35 additional cycles for 15s at 94°C, annealing for 30s at 72°C. A melting procedure with a heating rate of 0.5°C /10s was performed to create melting curves ranging from 95 °C. The gene expression levels were estimated using the 2 -ΔΔCt . Additional le 7 lists the sequences used for the primers.

Measuring the m 6 A Content
The overall content of mRNA m 6 A was measured by a methylation quanti cation kit of EpiQuik RNA (Epigentek, P-9005, NY, USA). In short, a standard curve was constructed at concentrations of 0.01-0.5 ng/ μL by positive control. The equivalent RNA solution (1-8 μL) and negative control were applied to the strip wells. The plate was wrapped with para lm, hatching for 1.5 h at 37 °C. Then, the wells were washed three times and added the 1:1000 diluted capture antibody at room temp for 1 h. After washing thrice, the detection antibody (1:2000 dilution) and enhancer solution were applied to every well hatched at room temp for 30 min. After ve washes, detection solutions were placed on each well and hatched for 10 minutes at room temperature to protect from light. Finally, stop solution applied to each wells and absorbance read with microplate reader at 450 nm.
MeRIP-seq and mRNA sequencing According to manufacturer's protocol, the total RNA was extracted using Trizol reagent (Invitrogen, CA, USA). Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 7.0 were used to evaluate the total RNA quality and quantity. Nearly over 200 ug total RNA was performed to isolate Poly (A) mRNA through magnetic beads (Invitrogen) attached poly-T oligo. After purifying, poly (A) mRNA fractions are broken into 100nt long oligonucleotides using divalent cations at high temperatures.
Then, the fragmentation of broken RNA were incubated in an IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5 % Igepal CA-630) supplied with BSA (0.5-1 μg / μl) for 2 hours at 4 ° C with m 6 A-speci c antibody (No. 202003, Synaptic Systems, Germany). The blend was further hatched with protein-A beads and eluted with the extraction buffer (1×IP buffer and 6.7 mM m 6 A), following precipitated by 75% ethanol.
DUTP method was used to transform eluted m 6 A-containing fragments (IP) and undisposed input control fragments into nal cDNA library according to a strand-speci c library preparation. The paired-end libraries of insert sizes were 100 ± 50bp. Afterwards, we implemented the paired-end sequencing on an Illumina Novaseq TM 6000 system at LC-BIO Bio-tech ltd (Hangzhou, China), using a suggested protocol from the manufacturer.

Sequencing data analysis
First of all, in-house perl scripts and Cutadapt [62] were performed eliminate the reads containing contaminants of the adapter, bases of low quality and indeterminate. Meanwhile, the quality of the sequence was validated using fastp. The reads were mapped to the Bos_mutus genome (Version: BosGru_v2.0) by HISAT2 [63] with default parameters. Using R package exomePeak [64] identify the m 6 A peaks from mapped reads of IP and input libraries with bed or bam format so as to con gure for viewing on IGV software (http://www.igv.org/) or the UCSC genome browser. De novo and de ned motif were identi ed by MEME [65] and HOMER [66], accompanied by perl scripts in the house seeking the motif with respect to peak. Called peaks were annotated using ChIPseeker [67] by intersection with gene architecture. StringTie [68] calculated the expression level of all mRNAs from input libraries, in which normalized with FPKM (FPKM=[total exon fragments/mapped reads (millions) ]). The differentially expressed mRNAs were collected by R package edgeR [69] with the |log2 (fold change)| > 1 and p-value < 0.05. GO seq R package was performed to the Gene Ontology (GO, http://www.geneontology.org/) enrichment analysis for the differentially expressed genes [70]. Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/) database is a major resource for learning high-level functions and utilities of biological systems. The statistical enrichment tests for genes of differential expression in the KEGG pathways were used the KOBAS software [71].

Statistical analysis
The SPSS 22 software package was used to evaluate statistics. A one-way test of variance assessed the signi cance of the differences between all of the groups. Statistically signi cant was the degree of probability * p<0.05;** p<0.01. Values are shown as mean ± SEM. Declarations System (CARS-37) supplied the source of experimental animals. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Authors' contributions
The experiments were conceived and designed by Y.Z. and P.Y. The experiments were performed by Y.Z.
The experiments were assisted with J.P., X.W., X.G., M.C., P.B., X.D., C.L. The paper was written by Y.Z. and revised by Q. K. All authors have read and agreed to the published version of the manuscript.    The abundance of m6A on KLF9 and FOXO1 mRNA transcripts observed by m6A-seq in preadipocyte and adipocyte. (c) GO terms of genes with differential m6A peaks between preadipocyte and adipocyte. (d) The top20 markedly enriched pathways for the genes of differential peaks between preadipocyte and adipocyte.