Skip to main content
Download PDF

Profiling of N6-methyladenosine methylation in porcine longissimus dorsi muscle and unravelling the hub gene ADIPOQ promotes adipogenesis in an m6A-YTHDF1–dependent manner

Journal of Animal Science and Biotechnology volume 14, Article number: 50 (2023) Cite this article

Abstract

Background

Intramuscular fat (IMF) content is a critical indicator of pork quality, and abnormal IMF is also relevant to human disease as well as aging. Although N6-methyladenosine (m6A) RNA modification was recently found to regulate adipogenesis in porcine intramuscular fat, however, the underlying molecular mechanisms was still unclear.

Results

In this work, we collected 20 longissimus dorsi muscle samples with high (average 3.95%) or low IMF content (average 1.22%) from a unique heterogenous swine population for m6A sequencing (m6A-seq). We discovered 70 genes show both differential RNA expression and m6A modification from high and low IMF group, including ADIPOQ and SFRP1, two hub genes inferred through gene co-expression analysis. Particularly, we observed ADIPOQ, which contains three m6A modification sites within 3′ untranslated and protein coding region, could promote porcine intramuscular preadipocyte differentiation in an m6A-dependent manner. Furthermore, we found the YT521‑B homology domain family protein 1 (YTHDF1) could target and promote ADIPOQ mRNA translation.

Conclusions

Our study provided a comprehensive profiling of m6A methylation in porcine longissimus dorsi muscle and characterized the involvement of m6A epigenetic modification in the regulation of ADIPOQ mRNA on IMF deposition through an m6A-YTHDF1-dependent manner.

Background

Intramuscular fat (IMF) content is a critical indicator in pork consume, and also linked to insulin resistance [1], aging [2] and obesity [3] in human. Pig works as an ideal human biomedical model with advantage over primates and other livestock due to its high similarities with human being from the anatomy and physiology [4]. Therefore, illustrating the molecular mechanism underlying the IMF deposition is vital for pork consumption and human health.

N6-Methyladenosine (m6A) is the most prevalent post-transcriptionally modification in eukaryotic cells, emerging as an important epigenetic regulator in various physiological processes [5, 6]. Dynamic mRNA m6A modification is regulated by dedicated methyltransferases (“writers”) and demethylases (“erasers”) [7]. RNA-binding proteins (“readers”) could recognize m6A-containing transcripts to drive RNA processes [8, 9], such as mRNA stability [9], splicing [10] or translation [11]. For instance, YT521‑B homology domain family protein 1 (YTHDF1) promotes breast cancer metastasis via enhancing FOXM1 translation in an m6A-dependent manner [12]. Fat mass and obesity-associated (FTO) protein regulates the splicing of adipogenic regulatory factor RUNX1T1 through affecting m6A level around splice site [13]. It has been reported that m6A is highly enriched around the stop codons and 3’UTRs [5]. Recent progress also indicated that m6A methylation of the 3’UTR of FLC causing depletion of its mRNA, controlling flowering in Arabidopsis [14].

Accumulating evidences suggested that m6A modification played important roles in regulating various aspects of mRNA metabolism during adipose tissue expansion [15,16,17,18]. For instance, NADP modulates m6A methylation and adipogenesis by enhancing FTO activity in 3 T3-L1 preadipocytes [19]. Consistently, Zfp217 mediates mRNA m6A methylation through FTO and YTHDF2 to regulate adipogenesis [20]. Furthermore, m6A modification of two adipogenesis-related genes, UCP2 and PNPLA2, would both regulate adipogenesis between Chinese indigenous breed Jinhua (fatty) and Western commercial breed Landrace (lean) in backfat, whereas in an opposite way [21]. Although it has been reported that YTHDF1 directly targets MTCH2 to promote adipogenesis in porcine intramuscular preadipocytes, our understanding about the function of m6A modification in IMF deposition was still limited.

Here we aimed to provide a valuable resource to determine the effects of m6A modified genes potentially involving in adipogenesis of IMF, permitting us to better understanding how to improve pork quality and providing potential target for therapy of obesity.

Materials and methods

Animal, phenotype and sample collection

This study utilized a mosaic swine population to uncover the relationship of m6A regulation mechanism and IMF deposition. The heterogeneous pig stock was derived from eight founder breeds (F0) consisting of the four Western commercial breeds (Duroc, Large White, Landrace and Pietrain pigs) and the four Chinese indigenous breeds (Erhualian, Laiwu, Bamaxiang and Tibetan pigs). All the pigs were raised under the same condition and purposeful mating, crossbreed strategy in detail was described previously [22, 23]. Animals were slaughtered in commer abattoir at 240 ± 10 d. We selected the longissimus dorsi muscle (LDM) from the 6th generation (F6; average IMF: 2.28%, range 0.92%–7.45%) [23]. LDM was obtained between the 3rd and 4th lumbar vertebrae, and flash frozen in liquid nitrogen and stored at −80 °C before use. The intramuscular fat content was measured using the routine Soxhlet extraction method [24].

Intramuscular preadipocytes cells were isolated from the LDM of 3-day-old Duroc-Landrace-Yorkshire piglets under sterile conditions [15]. The experimental procedures were in compliance with guidelines of the Committee on Animal Care and Use and Committee on the Ethic of Animal Experiments of Zhejiang University (Hangzhou, China).

RNA extraction and m6A RNA immunoprecipitation sequencing

Total RNA was isolated and purified using Trizol reagent (Invitrogen, Carlsbad, CA, USA) refer to the instruction, criteria with RIN > 7.0, total RNA > 50 μg, concentration > 50 ng/μL and OD260/280 > 1.8 were left for subsequent use. Poly (A) RNA is purified from 50 μg total RNA using DynabeadsTM Oligo (dT)25–61005 (Thermo Fisher Scientific Baltics UAB; Vilnius, Lithuania) using two rounds of purification. Then the poly(A) RNA was fragmented into small pieces using Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA) under 86 °C for 7 min.

Approximately 50 ng of fragmented mRNA was saved as input sample, which was used to eliminate the background. m6A-sepecific methylated RNA sequencing was performed according to the previous report [25]. In brief, the other fragmented mRNA was incubated with 3 μg methylated RNA-specific antibodies (No. 202003, Synaptic Systems, Göttingen, Germany) in RIP buffer (150 mmol/L NaCl, 10 mmol/L Tris and 0.1% NP-40) at 4 °C. After 2 h, adding the washed protein A/G magnetic beads (Millipore, Billerica, MA, USA) and incubating at 4 °C for further 2 h. Beads, washed 6 times in RIP buffer, incubated with immunoprecipitation buffer (Sigma-Aldrich, St Louis, MO, USA) to elute RNA. Immunoprecipitated RNA was extracted with phenol/chloroform, and RNA samples were sent for next-generation sequencing. All libraries were sequenced for 150 bp paired-end sequencing under an Illumina Novaseq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor’s recommended protocol.

Quantitative of m6A level by liquid chromatography-tandem mass spectrometry (LC-MS/MS)

Quantification of m6A in mRNA was conducted based on the previous study [26]. In brief, 300 ng of mRNA was digested by nuclease P1 (2 U) at 42 °C for 2 h, followed by the addition of alkaline phosphatase (0.5 U) with incubation at 37 °C for 2 h. The total amount of m6A in RNA was measured using Waters Acquity UPLC coupled to a Waters Xevo TQ mass spectrometer (Waters, Milford, USA). Quantification was achieved by comparing with the standard curve obtained from pure nucleoside standards. The ratio of m6A to A was calculated based on the determined concentration.

RNA mapping and quality control

Raw data were evaluated with FastQC v0.11.9 [27], the heading 10 bp were removed using trimmomatic v0.39 on account of GC bias [28]. Clean data were mapped to Sus scrofa 11.1 using STAR v2.7.8a, SAMtools v1.11 was used for sorting and marking duplicated reads [29, 30]. IP data were performed the same mapping procedure as input data.

MeRIP-seq data analysis

For IP data, m6A peak calling was conducted by MACS2 with “--nomodel -g 2.5e9 --broad --keep-dup all” on whole transcript level. Differentially peaks were identified with in-house R script according to previous study [31, 32]. Briefly, bedtools was used to combine all peaks from High and Low group into a reference peak. Normalized depth of each peak was inferred by following method: Normarlized depth = ((IP reads of Peak Region/Total reads of IP sample) − (Input reads of Peak Region/Total reads of input sample))/Length of peak. Total number of each sample’ read was calculated by SAMtools v1.11 flagstats based on BAM file. Coverage of peaks were inferred using SAMtools v1.11 bedcov.

Differentially methylated peaks (P < 0.05 and abs (log2foldchange ) > log21.5) were identified by comparing average normalized depth of each peak between High and Low group using t-test in R program. VEP software was using for annotating the differential peaks. HOMER software was applying for uncovering the motif in conserved peak regions.

Input data analysis

Input data were used for annotating, merging and quantifying with StringTie v2.1.7. raw counts of transcripts then were normalized (described in the legend of Fig. 3f) and low expression genes (gene counts > 9 in less than 4 samples) were filtered. Differential expression transcripts/genes were uncovered by the DESeq2 software [33, 34]. (abs (log2foldchange )) > log21.5 and P < 0.05 were identified as differentially expressing transcripts/genes.

Principal component analysis was conducted with DESeq2 [33]. Briefly, high expression gene counts were used for constructing DESeq data with the function DESeqDataSetFromMatrix(). And then data was normalized by the function rlogTransformation(). PCA was inferred with the function plotPCA() and visualized in R program. Pheatmap package was performed for visualing heatmap.

Weighted gene co-expression network analysis

Co-expression network analysis was performed with WGCNA (Wegithed Gene Co-expression Network Analysis) R package [35]. Briefly, raw count of genes were infered from the input data, and then low expression genes (gene counts > 9 in less than 4 samples) were filtered. Reserved gene counts were normalizd with transcript per millon (TPM) method. The soft threshold power β is determined based on the standard scale-free network, inferred from the function pickSoftThreshold().The adjacency matrix was calculated using topological overlap measure (TOM) [36], hierarchically clustering coexpressed genes into modules. Module-trait associations were calculated as the Pearson’s correlation between the module eigengene and trait of interest [37]. The most relevant traits of module was selected for analyzing their biological function and uncovering hub genes. Hub genes are a group of genes with the highest connectivity, and determine the characteristics of the gene module. We defined hub genes which are the significant correlation with clinical characteristics (Gene Significance, GS > 0.2) and high module characterization (Module Membership, MM > 0.8) in the module.

Functional enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted by ClueGO in Cytoscape v3.9.0. Pathways with P ≤ 0.05 were selected, P-value was chosen from the term P-value corrected with Bonferroni step down. GO ontologies involve biological process, cellular component and molecular function.

Western blot analysis

Cells were lysed with the mixture containing cell lysis buffer for Western and IP and 1% phenylmethanesulfonyl fluoride (PMSF) (Biosharp, Beijing, China) on ice to extract protein. Protein samples were separated by SDS-PAGE and then transferred to polyvinylidene difluoride membranes. And the membranes were blocked with 5% non-fat milk at room temperature for 1 h, then incubated with the primary antibody overnight at 4 °C and next with the secondary antibody for 1 h at room temperature. The protein bands were visualized using ECL Protect from Light (Biosharp) and quantified using Image J software. The primary antibodies used in this study were as follows: ADIPOQ (sc-136131, Santa Cruz, Watsonville, CA, USA, diluted 1:200), FLAG (20543–1-AP, Proteintech, Rosemont, IL, USA, diluted 1:1,000), YTHDF1 (17479–1-AP, Proteintech, diluted 1:1,000), GFP (ET160–25, Huabio, Hangzhou, China, diluted 1:5,000), β-actin (M1210–2, Huabio, diluted 1:5,000). The secondary antibodies were as follows: goat anti-mouse IgG-HRP (HA1006, Huabio, diluted 1:2,000), goat anti-rabbit IgG-HRP (HA1001, Huabio, diluted 1:2,000).

Real-time quantitative PCR (qPCR) analysis

Total RNA was extracted using TRIzol (Biosharp) according to the product protocol. After examination of RNA purity and concentration, 2 μg RNA was used as a template to reverse transcribe to cDNA by using M-MLV Reverse Transcriptase Kit (Invitrogen). Reverse transcription conditions were under 5 min at 25 °C, 45 min at 50 °C, 5 min at 85 °C. qPCR analysis was performed using the SYBR Green PCR Master Mix (Roche, Basel, Switzerland) with the ABI Step-One PlusTM Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). Relative level of RNA expression was determined with 2−ΔΔCt method after normalization to GAPDH. Reaction conditions were 95 °C for 1 min, 40 cycles of 95 °C for 15 s and 60 °C for 30 s. Primers used in this study were listed in Table 1.

Table 1 Primer sequences used in this work
Full size table

Oil Red O staining

Oil Red O staining was performed as following procedures: cells were washed and fixed with 10% formalin for 1 h, and then washed 3 times with 60% isopropanol. Cells were stained with Oil Red O working solution (0.35% Oil Red O dye in 60% isopropanol) for 10 min, and further washed 4 times with distilled water. Cells were eluted the stained lipid droplets using 100% isopropanol for 10 min, and then measuring optical density (OD) at 500 nm to conduct the quantitative of lipid content.

Intramuscular preadipocytes (IMF cells) isolation

IMF cells were isolated based on the previous study [15]. Briefly, the LDM of 3-day-old Duroc-Landrace-Yorkshire piglets were separated under sterile conditions. Visible connective tissue was removed and finely minced. Muscle tissues were then digested in a digestion buffer consisting of 1 mg/mL collagenase type I (Gibco, Carlsbad, CA, USA) in a shaking water bath for 1.5 h at 37 °C. The digested sample was filtered aseptically through 80 and 200 μm nylon mesh filters to isolate cells. Filtered cells were then washed 3 times with Dulbecco’s Modified Eagle Medium (DMEM) via centrifugation at 1,500 r/min for 5 min. Cells were seeded in growth medium that consisted of DMEM medium containing 10% fetal bovine serum (Gibco) and 1% penicillin-streptomycin (Gibco). After 1 h, cells were rinsed with DMEM medium to remove unadhered cells, and the adhered cells consisted of pure IMF cells.

Cell culture and adipocyte differentiation

Cells were cultured in DMEM containing 10% fetal bovine serum (Gibco) and 1% penicillin-streptomycin (Gibco). At 2 d after confluence, defined as d 0, cells were induced to differentiation medium containing 0.5 mmol/L 3-Isobutyl-1-methylxanthine (IBMX), 1 μmol/L dexamethasone and 5 μg/mL insulin (Sigma, St. Louis, MO, USA). On d 2, the medium was replaced with maintenance medium containing 5 μg/mL insulin (Sigma) every 2 d until d 8. Two hundred and ninety-three T cells were cultured in DMEM/F12 medium containing 10% fetal bovine serum and 1% penicillin-streptomycin (Gibco). Cells were uniformly cultured in a 5% CO2 incubator with 37 °C.

Cell transfection, plasmids and RNA knockdown

The plasmids and siRNA transfections were performed using Hieff Trans™ Liposomal Transfection Reagent and Hieff Trans™ in vitro siRNA/miRNA Transfection Reagent (Yeasen, Shanghai, China), according to the product protocol. The adenoviruses ADV4-ADIPOQ-CDS wild-type (ADV4-ADIPOQ-CDS-WT), ADV4-ADIPOQ-CDS mutant (m6A C534 and C570 were replaced by T, ADV4-ADIPOQ-CDS-MUT) and ADV4-ADIPOQ-CDS negative control (ADV4-ADIPOQ-CDS-NC) were generated by GenePharma (Shanghai, China). IMF cells were infected with the multiplicity of infection (MOI) of 25:1 by ADV4-ADIPOQ-CDS-WT, ADV4-ADIPOQ-CDS-MUT and ADV4-ADIPOQ-CDS-NC, respectively, and added 1 μg/mL polybrene to improve the infection efficiency, according to GenePharma’s protocol. Porcine YTHDF1 cDNA was generated via PCR and cloned into the pFLAG-CMV2 expression plasmid. Sequences of siRNA, synthesizd by GenePharma (Shanghai, China), were as follows: siADIPOQ-F, 5′- AGAAAGCGCCUAUGUCUACTT-3′ and siADIPOQ-R, 5′-GUAGACAUAGGCGCUUUCUCC-3′; siYTHDF1-F, 5′-UUAGUAUCCUGUCCUUUUGUU-3′ and siYTHDF1-R, 5′-CAAAAGGACAGGAUACUAAAG-3′.

m6A-specific methylated RNA immunoprecipitation real-time PCR

m6A-qPCR analysis was conducted according to previously report [38]. Briefly, mRNAs fragmented by RNA fragmentation reagent (Invitrogen) at 70 °C for 15 min. 10% of fragmented RNAs was used as input control mRNAs. The remaining 90% was immunoprecipitated with anti-m6A antibody coupled to Dynabeads (Invitrogen) in immunoprecipitation buffer (RNase inhibitor, 10 mmol/L Tris-HCl, 150 mmol/L NaCl, 0.1% Igepal CA-630 [Sigma]) at 4 °C for 2 h. mRNAs containing m6A were eluted twice with m6A 5′-monophosphate sodium salt (Sigma) at 4 °C for 1 h. After ethanol precipitation, all mRNAs were reversely transcribed into cDNA by M-MLV reverse transcriptase (Invitrogen). And then m6A enrichment was determined by qPCR. Data were analyzed with the 2−ΔΔCt method, and the relative enrichment of m6A in each sample was calculated by normalizing to input. The primers were as follows: ADIPOQ-CDS-F, 5′- TCCTTCCACATCACGGTCTACT-3′ and ADIPOQ-CDS-R, 5′- CTCCAGATAGAGGAGCACAGAG-3′; ADIPOQ-3’UTR-F, 5′-CCACTGTGTTTCCTCAGGTTC-3′ and ADIPOQ-3’UTR-R, 5′- CCACAGCCCTGTGTTTGACTT-3′.

RNA immunoprecipitation assay

The experiment pipeline was performed according to the previous research [39]. Briefly, FLAG-YTHDF1 overexpressed IMF cells were lysed in lysis buffer for 30 min at 4 °C and the supernatant was collected for further use. We saved 50-μL aliquot of cell lysate as input, and the remaining was incubated with anti-FLAG or immunoglobulin G (IgG) antibody-conjugated magnetic beads (Sigma) for 4 h at 4 °C. The beads were washed with buffer containing 0.1% SDS and proteinase K (Invitrogen), detecting fold enrichment with qPCR.

Dual-luciferase reporter and mutagenesis assays

To evaluate the effect of 3’UTR m6A site on ADIPOQ expression, wild type or mutant (m6A A650 was replaced by T) of ADIPOQ-3’UTR was inserted into downstream of pmirGLO Dual-Luciferase vector (Promega, Madison, WI, USA). After 48 h post transfection, the activities of firefly luciferase and Renilla luciferase in each 24-well plates’ well were determined by a Dual-Luciferase Reporter Gene Assay Kit (Yeasen) according to the product protocol.

Statistical analysis

All data were presented as mean ± SEM. Statistical differences in the dual luciferase reporter assay were determined by Mann-Whitney test, and other statistical significance between multiple groups were determined by Student’s t-test with GraphPad Prism 9. P < 0.05 was considered exceeding the significant level.

Results

Description of m6A modification between high and low IMF content groups

To investigate the role of m6A modification on adipogenesis in LDM, we collected 20 extreme phenotypic samples of IMF content from the 6th generation individuals in a unique heterogeneous swine population, which exhibits a large variation of IMF content [22, 23]. The samples were divided into high and low group according to IMF content (High and Low), and LC-MS/MS was performed to evaluate the m6A modifications levels across the samples. We found that IMF content (left in Fig. 1a) and level of m6A modifications (right in Fig. 1a) displayed opposite trend across the group, while both of those were significantly divergent among High and Low (P < 0.01), in agreement with previous study [15]. We uncovered 20,738 and 20,117 peaks among High and Low (Fig. 1b), respectively. A total of 23,250 peaks as a m6A modification panel within this population were obtained by “bedtools merge -d 0”. Conserved m6A modification motif among the panel was concordance with previous study (RRACH (R = G or A and H = A, C or U)) using HOMER (Fig. 1d). m6A modification sites were accumulating at the stop codon site (Fig. 1c) [15]. Peaks, annotated with ChIPseeker, were mainly enriched in the 3’UTR (Fig. 1e). These results together suggested our data was credible to further investigate the effect of m6A modification on lipid deposition in LDM.

Fig. 1
figure 1

Overview of m6A modification in High and Low IMF content groups. a Intramuscular fat ratio (right) and m6A/A content (left) among High and Low group, n = 10. b Venn diagram of peaks among two groups. c Density of m6A modification across mRNA region. d Conserved motif in m6A peaks using HOMER software. e Annotation of location of m6A peak at whole-transcript level. ***P < 0.001

Full size image

Identifying co-expression gene module of LDM

Co-expression network analysis enable us to identify genes which have a tendency to show a coordinated expression pattern among samples, uncovering the complexity of a cellular transcription network [37, 40]. Thus, we conducted the WGCNA software [35] to construct a co-expression network with 13,245 highly expressing genes (≥ 10 reads in at least 16 individuals) among 19 samples from the input RNA sequencing (RNA-seq) data of m6A-seq. L96 was excluded for outlier clustering according to the PCA and heatmap (Fig. S1a–d). We then chosen the optimal weighting coefficient β = 7 to construct the network based on pickSoftThreshold parameter in WGCNA. Figure 2a shows the cluster tree of the 19 samples and the corresponding traits information. Of 33 identified gene modules (Fig. 2b), MEdarkturquoise (module eigengene in dark turquoise) with 70 genes (Fig. 2c; Table S1) was detected significantly positively related to IMF content (r = 0.62; P = 0.004) and highly negatively associated with m6A content (r = − 0.51; P = 0.03) respectively, suggesting these genes within the module potentially participant in fat deposition. To investigate the underlying role of these co-expression genes, we performed the KEGG and GO pathway enrichment analysis by ClueGO in Cytoscape v3.9.0 [41]. The significant biological processes were involved in several adipogenesis related pathways, such as regulation of fat cell differentiation (ADIPOQ, BMP2, CEBPα, PPARγ, SFRP1), positive regulation of fat cell differentiation (BMP2, CEBPα, PPARγ, SFRP1) and PPAR signaling pathway (ADIPOQ, PLIN1, PPARγ) (Table S5). In addition, we identified 12 hub genes (ADIPOQ, PLIN1, UNC93A, SFRP1, HACD2, SNCG, SDR16C5, PPARγ, ITIH3, FFAR4, SORL1 and ACE2) from the dark turquoise module based on |geneModuleMembership| > 0.8 and |geneTraitSignificance| > 0.2 (Table 2). Of these, ADIPOQ, PLIN1, SFRP1, PPARγ and FFAR4 have been reported to participate in adipogenesis related function [42,43,44,45]. Remarkably, ADIPOQ, PLIN1 and FFAR4 were identified higher expression in subcutaneous fat and intramuscular fat compared with LDM among the same heterogeneous swine population [23], hinting these hub genes may play critical roles in adipogenesis.

Fig. 2
figure 2

Network analysis of MeRIP-seq input data of LDM samples. a Sample dendrogram from 19 LDM samples and trait heatmap including content of IMF and m6A modifications. Color intensity is directly proportional to the value of corresponding trait. b Cluster dendrogram of 13,245 highly expressing genes. Thirty-three co-expression modules were identified, each color represents a module. c Heatmap of the correlation between module eigengenes (MEs) and traits. Left value is correlation, and right enclosed in bracket is P-value. d Scatter plot for 70 genes in dark turquoise module, gene significance (GS) > 0.2 and module membership (MM) > 0.8 were selected as hub gene

Full size image
Table 2 Hub genes screened with WGCNA in porcine LDM
Full size table

ADIPOQ gene display significantly difference in both m6A modification and RNA expression between high and low group

To determine the role of m6A modification in intramuscular fat, we annotated the differential peak regions: 953 and 654 genes (Fig. 3a) were uniquely modified with m6A across the High and Low, respectively. One thousand and eighty-five genes (Fig. 3b; Table S2) were identified for significantly differential modified (abs (log2foldchange) > log21.5; P < 0.05). Gene ontology analysis of these m6A modified regions were significantly enriched in lipoprotein related functions (Fig. 3c), suggesting mRNA m6A in longissimus muscle play a potential role in regulating fat deposition. Among the 8 top significant differentially modified genes (according to P-value), we observed ADIOPQ (P = 5.09E−05) and SH3PXD2B (P = 1.28E−04) were reported to regulate the fat cell differentiation (Fig. 3b) [46, 47]. Similarly, we discovered 422 differential expression genes (abs (log2foldchange) > log21.5; P < 0.05) among the High and Low based on input RNA-seq data from m6A-seq (Fig. S1e; Table S3). Gene enrichment analysis revealed lipid droplet (CIDEC, PLIN1, PNPLA3, SDR16C5, TMEM135) and PPAR signaling pathway (ADIPOQ, AQP7, aP2, PLIN1, PPARγ) enriched in up regulation gene set (Fig. S1f; Table S5). We also observed the ADIPOQ gene displaying significantly differential RNA expression (P = 7.65E−14) among High and Low. Accumulating evidence indicated that mRNA m6A modification could mediate transcription regulation [12, 48]. Thus, to investigate whether m6A contributes to translation regulation in longissimus muscle, we overlapped the genes significantly difference both in the level of m6A modification and RNA expression between High and Low. Finally, we found 70 target co-differential genes (Fig. 3d, e), including ADIPOQ and SFRP1, which were related to several pathways such as PPAR signaling pathway (ADIPOQ, AQP7, aP2) (Table S4). SFRP1 gene has been reported that inhibits the Wnt/β-catenin signaling pathway, regulating the adipogenesis both in human and murine [49]. ADIPOQ gene is expressed specifically in adipose tissue [50], which exhibited higher expression in porcine fat tissues including subcutaneous fat and intramuscular fat than LDM in the same population [23]. To illustrate the mechanism of m6A modification on regulating the adipogenesis, we then chosen the hub gene ADIPOQ with remarkably methylated and RNA expression co-differential for further investigation (Fig. 3f, g).

Fig. 3
figure 3

RNA expression differentially and m6A modification differentially genes between High and Low. a Venn diagram of m6A modified genes across High and Low groups. b Volcano plot of m6A modified differential gene, P < 0.05 and fold change > 1.5 were marked as differentially methylated genes (blue and red), fold change value is calculated by High/Low. c GO and KEGG pathways of down (blue) and up (red) regulated m6A modification genes. d Venn diagram and (e) four quadrant diagram of methylated and RNA expression differential genes (P < 0.05 and fold change > 1.5) between High and Low group, 70 genes were observed significantly co-differential in e. f and g m6A methylation and mRNA expression of ADIPOQ gene between High and Low group, n = 8. Normalized read count was employed for comparing the level of m6A methylation between High and Low. Normalized read Count = SRN/ITR, SRN is site of reads number, while ITR is individual of total reads. SRN was counted using SAMtools v1.11 bedcov, ITR was inferred using SAMtools v1.11 flagstats based on BAM file. Input and IP data were both under the same pipeline of normalization. h Protein level of ADIPOQ between High and Low, n = 3

Full size image

ADIPOQ promotes adipogenesis of preadipocytes in vitro

To re-validate whether ADIPOQ gene regulates adipogenesis, intramuscular preadipocytes were isolated for adipogenic differentiation by the standard IBMX, dexamethasone, and insulin (MDI) cocktail (Fig. 4a). The lipid accumulation and mRNA expression levels of adipogenic genes (PPARγ, CEBPβ and aP2) were significantly increased after MDI induction (Fig. S2a, b). Simultaneously, the expression of ADIPOQ mRNA and protein were significantly increased from d 0 to 8 (Fig. S2c; Fig. 4b).

Fig. 4
figure 4

ADIPOQ promote adipogenesis of preadipocytes in vitro. a Workflow of porcine intramuscular adipocytes inducing in vitro. b Protein levels of ADIPOQ in intramuscular preadipocytes at 0, 2, 4 and 8 d during adipogenesis. c The mRNA levels of ADIPOQ (48 h) after siRNA transfection of porcine intramuscular preadipocytes, n = 3. d The protein expression levels of ADIPOQ (48 h) after siRNA transfection of porcine intramuscular preadipocytes. e–g TAG content and Oil Red O staining of siADIPOQ at 8 d after adipogenic induction, n = 3. h RT-qPCR of PPARγ, CEBPβ and aP2 of siADIPOQ at 8 d after adipogenic induction, n = 3. i The mRNA expression levels of ADIPOQ after overexpression ADIPOQ (48 h), n = 3. j The protein expression levels of ADIPOQ after overexpression ADIPOQ (48 h). k, l TAG content and Oil Red O staining of ADIPOQ-overexpression at 8 d after adipogenic induction, n = 3. m RT-qPCR of PPARγ, CEBPβ and aP2 of ADIPOQ overexpression at 8 d after adipogenic induction, n = 3. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Full size image

Previous research had indicated that interference with ADIPOQ gene expression could inhibit the differentiation of porcine preadipocytes [42]. Thus, we established siRNA and overexpression plasmid to address the function of ADIPOQ in the process of adipogenic differentiation in our work. Not surprisingly, mRNA expression and protein level of ADIPOQ were significantly inhibited after siRNA interference at d 8 (Fig. 4c, d). Meanwhile, lipid accumulation of siADIPOQ was remarkably decreased according to triacylglycerol (TAG), Oil Red O staining and adipogenic genes (PPARγ, CEBPβ and aP2) mRNA expression (Fig. 4e–h). Next, we observed the overexpression of ADIPOQ in porcine intramuscular preadipocytes cell could significantly increase the levels of its mRNA expression and protein (Fig. 4i, j), promoting lipid accumulation (Fig. 4k–m). On the basis of above results, we concluded that ADIPOQ was expressed at the later stage of induction and promoted porcine intramuscular preadipocyte differentiation and lipid accumulation.

mRNA m6A modification can promote ADIPOQ expression

Although we acquired that ADIPOQ could promote the lipid accumulation in intramuscular preadipocytes, the role of m6A modification in ADIPOQ remain unclear [42, 47]. To further explore the function of mRNA m6A modification on ADIPOQ expression, we firstly scanned the transcript to uncover the m6A sites of ADIPOQ gene based on the RRACH conversed feature. Three potential m6A sites including one in 3’UTR (AGACT, chr13:124,645,333–124,645,337) and two in CDS (GGACA, chr13:124,644,484–124,644,488; GGACA chr13:124,644,520–124,644,524) were found in the longest ADIPOQ transcript ENSSSCT00000047495 (Fig. 3f). To explore the role of m6A modification in 3’UTR and CDS of ADIPOQ, we constructed the dual-luciferase reporter plasmid and adenovirus vector with mutation in 3’UTR (ADIPOQ-3’UTR-MUT) and CDS (ADIPOQ-CDS-MUT), respectively (Fig. 5a, b; Table S6). Analysis of m6A-IP-qPCR found that m6A methylation levels of ADIPOQ-CDS-WT and ADIPOQ-3’UTR-WT were higher than ADIPOQ-CDS-MUT and ADIPOQ-3’UTR-MUT, respectively (Fig. 5a). Luciferase assays results indicated that mutation of ADIPOQ 3’UTR significantly decreased the luciferase activity in 293 T cells (Fig. 5c). Consistently, the mRNA expression and protein level of ADIPOQ in ADIPOQ-CDS-WT IMF cells were also higher than ADIPOQ-CDS-MUT (Fig. 5d, e). We also found ADIPOQ-CDS-MUT decreases lipid accumulation (Fig. 5f, g) and adipocyte differentiation-related gene expression including PPARγ, CEBPβ and aP2, relative to ADIPOQ-CDS-WT (Fig. 5h). Taken together, we concluded that the m6A modification of ADIPOQ in 3’UTR and CDS could both promote its expression.

Fig. 5
figure 5

ADIPOQ promotes adipogenesis of preadipocytes in a m6A-dependent manner. a m6A-IP-qPCR analysis of ADIPOQ-3’UTR WT or MUT (A to T mutation) in 293 T cells, n = 3. b m6A-IP-qPCR analysis ADIPOQ-CDS WT or MUT (C to T mutation) in porcine intramuscular preadipocytes, n = 3. c Relative luciferase activity of WT or MUT of ADIPOQ-3’UTR in 293 T cells, n = 3. d The mRNA expression levels of ADIPOQ of porcine intramuscular preadipocytes with NC, WT or MUT of ADIPOQ-CDS, n = 3. e The protein expression levels of ADIPOQ of porcine intramuscular preadipocytes with NC, WT or MUT of ADIPOQ-CDS. f and g TAG content and Oil Red O staining of porcine intramuscular preadipocytes with NC, WT or MUT of ADIPOQ-CDS at 8 d after adipogenic induction, n = 3. h RT-qPCR of PPARγ, CEBPβ and aP2 of porcine intramuscular preadipocytes with NC, WT or MUT of ADIPOQ-CDS at 8 d after adipogenic induction, n = 3. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Full size image

YTHDF1 mediates the regulation of ADIPOQ in an m6A-dependent manner

We then explored the mechanism about how m6A modification regulated ADIPOQ expression. YTHDF1 was reported to promote translation of m6A methylated transcripts [51]. Regarding the m6A sites in 3’UTR or CDS could promote the translation of ADIPOQ, we assumed that ADIPOQ is the target of YTHDF1. Thus, we performed YTHDF1 knockdown and overexpressing experiments to identify whether it could regulate ADIPOQ expression. Not surprisingly, YTHDF1 knockdown decreased ADIPOQ protein expression (Fig. 6a), while YTHDF1 overexpression increased ADIPOQ protein expression (Fig. 6b). RIP-qPCR assay revealed that ADIPOQ interacted with YTHDF1-FLAG, which confirmed that ADIPOQ is the target of YTHDF1 (Fig. 6c, d). To further explore whether YTHDF1 targets and recognizes the ADIPOQ mRNA m6A modification site, we transferred YTHDF1 overexpression plasmid into ADIPOQ-3’UTR-WT (or MUT) and ADIPOQ-CDS-WT (or MUT) cells, respectively. Overexpression of YTHDF1 increased luciferase activity and ADIPOQ protein level in ADIPOQ-3’UTR-WT 293 T cells, but no change in ADIPOQ-3’UTR-MUT (Fig. 6e, f) cells. Similarly, overexpressing YTHDF1 increased the mRNA and protein expression of ADIPOQ in ADIPOQ-CDS-WT IMF cells but no change in ADIPOQ-CDS-MUT cells (Fig. 6g, h). Moreover, we also observed overexpressing YTHDF1 increases lipid accumulation (Fig. 6i, j) and adipocyte differentiation-related gene expression including PPARγ, CEBPβ and aP2 (Fig. 6k) in ADIPOQ-CDS-WT but not in ADIPOQ-CDS-MUT. Collectively, these results together suggest YTHDF1 promotes the translation of hub gene ADIPOQ by recognizing m6A sites in both 3’UTR and CDS.

Fig. 6
figure 6

YTHDF1 regulate the translation of ADIPOQ in IMF cells. a The protein levels of ADIPOQ in porcine intramuscular preadipocytes transfected with siControl or siYTHDF1 (48 h). b The protein levels of ADIPOQ after overexpression YTHDF1 (48 h). c The protein levels of YTHDF1 of porcine intramuscular preadipocytes transfected with control or YTHDF1-FLAG plasmid (48 h). d RIP analysis of the interaction of ADIPOQ with FLAG in porcine intramuscular preadipocytes transfected with YTHDF1-FLAG plasmid. Enrichment of ADIPOQ with FLAG was measured by qPCR and normalized to input. e The protein levels of YTHDF1 of 293 T cells transfected with WT or MUT of ADIPOQ-3’UTR or YTHDF1 overexpression plasmid (48 h). f Relative luciferase activity of WT or MUT of ADIPOQ-3’UTR or YTHDF1 overexpression in 293 T cells, n = 3. g The mRNA expression levels of ADIPOQ of porcine intramuscular preadipocytes with WT or MUT of ADIPOQ-CDS or YTHDF1-overexpression plasmid (48 h), n = 3. h The protein expression levels of ADIPOQ of porcine intramuscular preadipocytes with WT or MUT of ADIPOQ-CDS or YTHDF1-overexpression plasmid (48 h). i and j TAG content and Oil Red O staining of porcine intramuscular preadipocytes with WT or MUT of ADIPOQ-CDS or YTHDF1 overexpression plasmid at 8 d after adipogenic induction, n = 3. k RT-qPCR of PPARγ, CEBPβ and aP2 of porcine intramuscular preadipocytes with WT or MUT of ADIPOQ-CDS or YTHDF1 overexpression plasmid at 8 d after adipogenic induction, n = 3. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Full size image

Discussion

In this work, we performed the m6A-seq of LDM from a unique heterogenous swine population to investigate the underlying mechanism of mRNA m6A modification regulating IMF deposition. We revealed that hub gene ADIPOQ could promote its mRNA translation in an m6A-YTHDF1-dependent manner, providing novel evidence of m6A methylation regulating adipogenesis.

Fat deposition is highly relevant to human health [52, 53], uncovering the mechanism porcine intramuscular adipogenesis is better for understanding gene regulation underlying the fat deposition of corresponding tissues in humans. Accumulating evidences demonstrate that m6A modification is involving in adipogensis pathway [17,18,19]. Although previous finding has been revealed that m6A modification of MTCH2 promotes adipogenesis in LDM when comparing obese Asian domesticated Jinhua pig and lean Western commercial pig, these results was still limited because of the selected validation gene merely obtain from top 10 methylation in Jinhua breed [15]. In this study, we possessed different hallmark from previous studies in that a unique swine population was used [15, 54], and found that some individuals exhibits a large variation of IMF content. More importantly, IMF content was negatively related to the mRNA m6A level across the High and Low group (P < 0.01), indicating a potential role of m6A in intramuscular fat deposition, which was consistent to previous study.

To explore the underlying role and mechanism of m6A in porcine intramuscular fat, we performed a large sample size of MeRIP data (n = 10 per group), which allowed us to explore more significant differential m6A modification sites. Tao et al. uncovered 5,872 and 2,826 m6A peaks respectively, in the porcine muscle and adipose tissue transcriptomes [54]. Here, we identified a total of 23,250 m6A peaks in this population, to our knowledge, it is largest m6A data set in procine intramuscular fat. Besides, the consensus motif sequence RRACH in our study was consistent with previous work. mRNA m6A sites were enriched around stop codons, sharing a smiliar distribution to those of human, mice and plants [55,56,57]. Taken together, using larger sample size and stringent m6A calling parameters, our results allow a reliable picture of the mRNA m6A epi-transcriptome in porcine skeletal muscles.

To uncover which key genes regulate adipogenesis in m6A-dependent manner, we performed gene co-expression network using WGCNA to explore the biologically relevant associations between phenotype and module [35]. Finally, we uncovering 70 modules among 19 high expression RNA input data, including 12 hub genes, were significantly corelated with IMF content and m6A modification level. Emerging evidences have indicated that WGCNA could reveal potential candidate gene in affecting the IMF content of Duroc [58] and Italian Large White pigs [59]. Thus, we overlapped RNA expression differential and m6A methylated differential genes, discovering 70 co-differential genes. We further found 2 hub gene ADIPOQ and SFRP1 including in co-differential gene set. These results largely advanced our knowledge towards co-expression networks in IMF deposition.

In this work, we found ADIPOQ gene displayed remarkably differential both in RNA expression (P = 7.65E−14) and m6A methylation (P = 5.09E−05). ADIPOQ has been identified as candidate gene for the metabolic syndrome and T2DM by genome wide associated study [60, 61]. Previous work also indicated ADIPOQ exhibited higher expression in both intramuscular fat and subcutaneous fat than LDM in the same swine population [23]. Consistently, previous study provided supportive evidence for silencing of ADIPOQ efficiently suppresses preadipocyte differentiation in porcine [42]. By establishing the lipogenesis model in vitro, we revalidated the ADIPOQ gene could promote the adipogenesis of porcine preadipocyte. We further found mRNA m6A modification could promote the expression of ADIPOQ and lipid accumulation by constructing the dual-luciferase reporter plasmid and adenovirus vector in 3’UTR and CDS, respectively.

Various m6A binding proteins, especially YTHDF family, have been proved their functions in different aspects, such as RNA translation, splicing, export or degradation [62, 63]. YTHDF1 selectively recognizes m6A in cytosolic mRNAs, recruiting initiation factor eIF3 to facilitate mRNA translation [51]. YTHDF2 brings m6A-modified translatable mRNAs to mRNA decay sites (e.g., P-bodies), and recruiting CC chemokine receptor 4-negative regulator of transcription complex to trigger mRNA deadenylation [9, 64]. YTHDF3 promotes mRNA translation in synergy with YTHDF1 and accelerated decay of m6A-containing mRNAs through interaction with YTHDF2. Accumulating evidences suggest YTHDF1 promote RNA expression via recognizing mRNA m6A site [65, 66]. YTHDF1 interacting with MTCH2 mRNA to enhance translation of its protein in porcine intramuscular preadipocytes [15, 67]. Thus, we here have conducted interference and overexpression YTHDF1 to confirm its function. Not surprising, we observed YTHDF1 promoting the translation of hub gene ADIPOQ, confirming that ADIPOQ was a target of YTHDF1 through m6A-IP and RIP experiments.

Conclusions

In conclusion, our study characterized the m6A modification genes which were potentially involved in regulating IMF deposition. Furthermore, we presented a novel regulatory mechanism of IMF deposition via the m6A-YTHDF1-ADIPOQ axis, highlighting the critical role of mRNA m6A modification of the hub gene in IMF adipogenesis.

Abbreviations

ACE2:

Angiotensin converting enzyme 2

ADIPOQ:

Adiponectin

aP2:

Fatty acid binding protein 4

AQP7:

Aquaporin 7

BAM:

Binary Alignment/Map format

BMP2:

Bone morphogenetic protein type 2

CDS:

Coding sequence

CEBPα:

CCAAT enhancer binding protein alpha

CEBPβ:

CCAAT enhancer binding protein beta

CIDEC:

Cell death inducing DFFA like effector c

DMEM:

Dulbecco’s modified eagle medium

FFAR4:

Free fatty acid receptor 4

GO:

Gene ontology

HACD2:

3-hydroxyacyl-CoA dehydratase 2

IBMX:

3-Isobutyl-1-methylxanthine

IgG:

Immunoglobulin G

IMF:

Intramuscular fat

IP:

Immunoprecipitation

ITIH3:

Inter-alpha-trypsin inhibitor heavy chain 3

KEGG:

Kyoto Encyclopedia of Genes and Genomes

LDM:

Longissimus dorsi muscle

m6A:

N6-methyladenosine

MDI:

Dexamethasone, and insulin

MTCH2:

Mitochondrial carrier homolog 2

MUT:

Mutation

OD:

Optical density

PLIN1:

Perilipin 1

PMSF:

Phenylmethanesulfonyl fluoride

PNPLA3:

Patatin like phospholipase domain containing 3

PPAR:

Peroxisome proliferator activated receptor

PPARγ:

Peroxisome proliferator activated receptor gamma

Qpcr:

Real-time quantitative PCR

RIP:

RNA immunoprecipitation

SDR16C5:

Short chain dehydrogenase/reductase family 16C member 5

SFRP:

Secreted frizzled-related protein 1

SH3PXD2B:

SH3 and PX domain-containing protein 2B

SNCG:

Synuclein gamma

SORL1:

Sortilin related receptor 1

TAG:

Triacylglycerol

T2DM:

Type 2 diabetes mellitus

TMEM135:

Transmembrane protein 135

UTR:

Untranslated region

UNC93A:

Unc-93 homolog A

WGCNA:

Weighted correlation network analysis

WT:

Wide type

YTHDF1–3:

YT521‑B homology domain family proteins 1–3

References

  1. Hilton TN, Tuttle LJ, Bohnert KL, Mueller MJ, Sinacore DR. Excessive adipose tissue infiltration in skeletal muscle in individuals with obesity, diabetes mellitus, and peripheral neuropathy: association with performance and function. J Physical therapy. 2008;88(11):1336–44. https://doi.org/10.2522/ptj.20080079.

    Article  Google Scholar 

  2. Schwenzer NF, Martirosian P, Machann J, Schraml C, Steidle G, Claussen CD, et al. Aging effects on human calf muscle properties assessed by MRI at 3 tesla. J Magn Reson Imaging. 2009;29(6):1346–54. https://doi.org/10.1002/jmri.21789.

    Article  PubMed  Google Scholar 

  3. Malenfant P, Joanisse D, Theriault R, Goodpaster B, Kelley D, Simoneau J. Fat content in individual muscle fibers of lean and obese subjects. Int J Obesity. 2001;25(9):1316–21. https://doi.org/10.1038/sj.ijo.0801733.

    Article  CAS  Google Scholar 

  4. Lunney JK, Van Goor A, Walker KE, Hailstock T, Franklin J, Dai C. Importance of the pig as a human biomedical model. Sci Transl Med. 2021;13(621):eabd5758. https://doi.org/10.1126/scitranslmed.abd5758.

    Article  CAS  PubMed  Google Scholar 

  5. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017;18(1):31–42. https://doi.org/10.1038/nrm.2016.132.

    Article  CAS  PubMed  Google Scholar 

  6. Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell. 2017;169(7):1187–200. https://doi.org/10.1016/j.cell.2017.05.045.

  7. Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m6A RNA methylation. Nat Rev Genet. 2014;15(5):293–306. https://doi.org/10.1038/nrg3724.

  8. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485(7397):201–6. https://doi.org/10.1038/nature11112.

  9. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N 6-methyladenosine-dependent regulation of messenger RNA stability. Nature. 2014;505(7481):117–20. https://doi.org/10.1038/nature12730.

    Article  CAS  PubMed  Google Scholar 

  10. Xiao W, Adhikari S, Dahal U, Chen Y-S, Hao Y-J, Sun B-F, et al. Nuclear m6A reader YTHDC1 regulates mRNA splicing. Mol Cell. 2016;61(4):507–19. https://doi.org/10.1016/j.molcel.2016.01.012.

  11. He PC, He C. m6A RNA methylation: from mechanisms to therapeutic potential. EMBO J. 2021;40(3):e105977. https://doi.org/10.15252/embj.2020105977.

  12. Chen H, Yu Y, Yang M, Huang H, Ma S, Hu J, et al. YTHDF1 promotes breast cancer progression by facilitating FOXM1 translation in an m6A-dependent manner. Cell Biosci. 2022;12(1):1–16. https://doi.org/10.1186/s13578-022-00759-w.

  13. Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, et al. FTO-dependent demethylation of N6-methyladenosine regulates mRNA splicing and is required for adipogenesis. Cell Res. 2014;24(12):1403–19. https://doi.org/10.1038/cr.2014.151.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Gong T, Han H, Tan Z, Ning Z, Qiao H, Yu M, et al. Segmentation and differentiation of periventricular and deep white matter hyperintensities in 2D T2-FLAIR MRI based on a cascade U-net. Front Neurol. 2022;13:1021477. https://doi.org/10.3389/fneur.2022.1021477.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Jiang Q, Sun B, Liu Q, Cai M, Wu R, Wang F, et al. MTCH2 promotes adipogenesis in intramuscular preadipocytes via an m6A-YTHDF1-dependent mechanism. FASEB J. 2019;33(2):2971–81. https://doi.org/10.1096/fj.201801393RRR.

  16. Wu R, Guo G, Bi Z, Liu Y, Zhao Y, Chen N, et al. m6A methylation modulates adipogenesis through JAK2-STAT3-C/EBPβ signaling. Biochim Biophys Acta Gene Regul Mech. 2019;1862(8):796–806. https://doi.org/10.1016/j.bbagrm.2019.06.008.

  17. Wang X, Wang Y. From histones to RNA: role of methylation in signal proteins involved in Adipogenesis. Curr Protein Pept Sci. 2017;18(6):589–98. https://doi.org/10.2174/1389203717666160627082444.

    Article  CAS  PubMed  Google Scholar 

  18. Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, et al. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther. 2021;6(1):74. https://doi.org/10.1038/s41392-020-00450-x.

  19. Wang L, Song C, Wang N, Li S, Liu Q, Sun Z, et al. NADP modulates RNA m(6)a methylation and adipogenesis via enhancing FTO activity. Nat Chem Biol. 2020;16(12):1394–402. https://doi.org/10.1038/s41589-020-0601-2.

  20. Song T, Yang Y, Wei H, Xie X, Lu J, Zeng Q, et al. Zfp217 mediates m6A mRNA methylation to orchestrate transcriptional and post-transcriptional regulation to promote adipogenic differentiation. Nucleic Acids Res. 2019;47(12):6130–44. https://doi.org/10.1093/nar/gkz312.

  21. Wang X, Sun B, Jiang Q, Wu R, Cai M, Yao Y, et al. mRNA m(6)a plays opposite role in regulating UCP2 and PNPLA2 protein expression in adipocytes. Int J Obes. 2018;42(11):1912–24. https://doi.org/10.1038/s41366-018-0027-z.

    Article  CAS  Google Scholar 

  22. Yang H, Wu J, Huang X, Zhou Y, Zhang Y, Liu M, et al. ABO genotype alters the gut microbiota by regulating GalNAc levels in pigs. Nature. 2022;606(7913):358–67. https://doi.org/10.1038/s41586-022-04769-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Zhang Y, Sun Y, Wu Z, Xiong X, Zhang J, Ma J, et al. Subcutaneous and intramuscular fat transcriptomes show large differences in network organization and associations with adipose traits in pigs. Sci China Life Sci. 2021;64(10):1732–46. https://doi.org/10.1007/s11427-020-1824-7.

    Article  CAS  PubMed  Google Scholar 

  24. Cameron N, Enser M, Nute G, Whittington F, Penman J, Fisken A, et al. Genotype with nutrition interaction on fatty acid composition of intramuscular fat and the relationship with flavour of pig meat. Meat Sci. 2000;55(2):187–95. https://doi.org/10.1016/s0309-1740(99)00142-4.

    Article  CAS  PubMed  Google Scholar 

  25. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell. 2012;149(7):1635–46. https://doi.org/10.1016/j.cell.2012.05.003.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Shafik AM, Zhang F, Guo Z, Dai Q, Pajdzik K, Li Y, et al. N6-methyladenosine dynamics in neurodevelopment and aging, and its potential role in Alzheimer's disease. Genome Biol. 2021;22(1):1–19. https://doi.org/10.1186/s13059-020-02249-z.

    Article  CAS  Google Scholar 

  27. Andrews S. FastQC: A quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

  28. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. https://doi.org/10.1093/bioinformatics/bts635.

    Article  CAS  PubMed  Google Scholar 

  30. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9. https://doi.org/10.1093/bioinformatics/btp352.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Zhu Y, Zhou Z, Huang T, Zhang Z, Li W, Ling Z, et al. Mapping and analysis of a spatiotemporal H3K27ac and gene expression spectrum in pigs. Sci China Life Sci. 2022;65(8):1517–34. https://doi.org/10.1007/s11427-021-2034-5.

    Article  CAS  PubMed  Google Scholar 

  32. Jiang T, Ling Z, Zhou Z, Chen X, Chen L, Liu S, et al. Construction of a transposase accessible chromatin landscape reveals chromatin state of repeat elements and potential causal variant for complex traits in pigs. J Anim Sci Biotechnol. 2022;13(1):112. https://doi.org/10.1186/s40104-022-00767-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21. https://doi.org/10.1186/s13059-014-0550-8.

    Article  CAS  Google Scholar 

  34. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5. https://doi.org/10.1038/nbt.3122.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. https://doi.org/10.1186/1471-2105-9-559.

  36. Yip AM, Horvath S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinformatics. 2007;8:22. https://doi.org/10.1186/1471-2105-8-22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Van Dam S, Vosa U, van der Graaf A, Franke L, de Magalhaes JP. Gene co-expression analysis for functional classification and gene–disease predictions. Brief Bioinform. 2018;19(4):575–92. https://doi.org/10.1093/bib/bbw139.

    Article  CAS  PubMed  Google Scholar 

  38. Dominissini D, Moshitch-Moshkovitz S, Amariglio N, Rechavi G. Transcriptome-wide mapping of N(6)-Methyladenosine by m(6)A-Seq. Methods Enzymol. 2015;560:131–47. https://doi.org/10.1016/bs.mie.2015.03.001.

    Article  CAS  PubMed  Google Scholar 

  39. Peritz T, Zeng F, Kannanayakal TJ, Kilk K, Eiriksdottir E, Langel U, et al. Immunoprecipitation of mRNA-protein complexes. Nat Protoc. 2006;1(2):577–80. https://doi.org/10.1038/nprot.2006.82.

    Article  CAS  PubMed  Google Scholar 

  40. Long J, Huang S, Bai Y, Mao J, Wang A, Lin Y, et al. Transcriptional landscape of cholangiocarcinoma revealed by weighted gene coexpression network analysis. Brief Bioinform. 2021;22(4):bbaa224. https://doi.org/10.1093/bib/bbaa224.

    Article  CAS  PubMed  Google Scholar 

  41. Gabriela B, Bernhard M, Hubert H, Pornpimol C, Marie T, Amos K, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):1091–3. https://doi.org/10.1093/bioinformatics/btp101.

    Article  CAS  Google Scholar 

  42. Gao Y, Li F, Zhang Y, Dai L, Jiang H, Liu H, et al. Silencing of ADIPOQ efficiently suppresses preadipocyte differentiation in porcine. Cell Physiol Biochem. 2013;31(2–3):452–61. https://doi.org/10.1159/000343381.

    Article  CAS  PubMed  Google Scholar 

  43. Sun Y, Zhai G, Li R, Zhou W, Li Y, Cao Z, et al. Rxrα positively regulates expression of the chicken plin1 gene in a pparγ-independent manner and promotes adipogenesis. Front Cell Dev Biol. 2020;8:349. https://doi.org/10.3389/fcell.2020.00349.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Aprile M, Ambrosio M, D'esposito V, Beguinot F, Formisano P, Costa V, et al. PPARG in human adipogenesis: differential contribution of canonical transcripts and dominant negative isoforms. PPAR Res. 2014;2014:537865. https://doi.org/10.1155/2014/537865.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Hilgendorf KI, Johnson CT, Mezger A, Rice SL, Norris AM, Demeter J, et al. Omega-3 fatty acids activate ciliary FFAR4 to control adipogenesis. Cell. 2019;179(6):1289–05. https://doi.org/10.1016/j.cell.2019.11.005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Fan C, Dong H, Yan K, Shen W, Wang C, Xia L, et al. Genome-wide screen of promoter methylation identifies novel markers in diet-induced obese mice. Nutr Hosp. 2014;30(1):42–52. https://doi.org/10.3305/nh.2014.30.1.7521.

    Article  CAS  PubMed  Google Scholar 

  47. Choi Y, Davis ME, Chung H. Effects of genetic variants in the promoter region of the bovine adiponectin (ADIPOQ) gene on marbling of Hanwoo beef cattle. Meat Sci. 2015;105:57–62. https://doi.org/10.1016/j.meatsci.2015.02.014.

    Article  CAS  PubMed  Google Scholar 

  48. Liu T, Wei Q, Jin J, Luo Q, Liu Y, Yang Y, et al. The m6A reader YTHDF1 promotes ovarian cancer progression via augmenting EIF3C translation. Nucleic Acids Res. 2020;48(7):3816–31. https://doi.org/10.1093/nar/gkaa048.

  49. Lagathu C, Christodoulides C, Tan CY, Virtue S, Laudes M, Campbell M, et al. Secreted frizzled-related protein 1 regulates adipose tissue expansion and is dysregulated in severe obesity. Int J Obes. 2010;34(12):1695–05. https://doi.org/10.1038/ijo.2010.107.

    Article  CAS  Google Scholar 

  50. Hu E, Liang P, Spiegelman BM. AdipoQ is a novel adipose-specific gene dysregulated in obesity (∗). J Biol Chem. 1996;271(18):10697–03. https://doi.org/10.1074/jbc.271.18.10697.

    Article  CAS  PubMed  Google Scholar 

  51. Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, et al. N(6)-methyladenosine modulates messenger RNA translation efficiency. Cell. 2015;161(6):1388–99. https://doi.org/10.1016/j.cell.2015.05.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Yang Y, Ding L, Zou X, Shen Y, Hu D, Hu X, et al. Visceral adiposity and high intramuscular fat deposition independently predict critical illness in patients with SARS-CoV-2. Obesity (Silver Spring). 2020;28(11):2040–8. https://doi.org/10.1002/oby.22971.

    Article  CAS  PubMed  Google Scholar 

  53. Cnop M, Landchild MJ, Vidal J, Havel PJ, Knowles NG, Carr DR, et al. The concurrent accumulation of intra-abdominal and subcutaneous fat explains the association between insulin resistance and plasma leptin concentrations : distinct metabolic effects of two fat compartments. Diabetes. 2002;51(4):1005–15. https://doi.org/10.2337/diabetes.51.4.1005.

    Article  CAS  PubMed  Google Scholar 

  54. Tao X, Chen J, Jiang Y, Wei Y, Chen Y, Xu H, et al. Transcriptome-wide N (6) -methyladenosine methylome profiling of porcine muscle and adipose tissues reveals a potential mechanism for transcriptional regulation and differential methylation pattern. BMC Genomics. 2017;18(1):336. https://doi.org/10.1186/s12864-017-3719-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Xiong X, Hou L, Park YP, Molinie B, Consortium GT, Gregory RI, et al. Genetic drivers of m(6)a methylation in human brain, lung, heart and muscle. Nat Genet. 2021;53(8):1156–65. https://doi.org/10.1038/s41588-021-00890-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Xu T, Wu X, Wong CE, Fan S, Zhang Y, Zhang S, et al. FIONA1-mediated m(6) a modification regulates the floral transition in Arabidopsis. Adv Sci (Weinh). 2022;9(6):e2103628. https://doi.org/10.1002/advs.202103628.

    Article  CAS  PubMed  Google Scholar 

  57. Tang Y, Chen K, Song B, Ma J, Wu X, Xu Q, et al. m6A-atlas: a comprehensive knowledgebase for unraveling the N6-methyladenosine (m6A) epitranscriptome. Nucleic Acids Res. 2021;49(D1):D134–43. https://doi.org/10.1093/nar/gkaa692.

  58. Zhao X, Hu H, Lin H, Wang C, Wang Y, Wang J. Muscle transcriptome analysis reveals potential candidate genes and pathways affecting intramuscular fat content in pigs. Front Genet. 2020;11:877. https://doi.org/10.3389/fgene.2020.00877.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Zappaterra M, Gioiosa S, Chillemi G, Zambonelli P, Davoli R. Muscle transcriptome analysis identifies genes involved in ciliogenesis and the molecular cascade associated with intramuscular fat content in large white heavy pigs. PLoS One. 2020;15(5):e0233372. https://doi.org/10.1371/journal.pone.0233372.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Siitonen N, Pulkkinen L, Lindström J, Kolehmainen M, Eriksson JG, Venojärvi M, et al. Association of ADIPOQ gene variants with body weight, type 2 diabetes and serum adiponectin concentrations: the Finnish diabetes prevention study. BMC Med Genet. 2011;12:5. https://doi.org/10.1186/1471-2350-12-5.

  61. Torres-Castillo N, Campos-Perez W, Rodriguez-Echevarria R, Rodriguez-Reyes SC, Martinez-Lopez E. A metabolically unhealthy phenotype is associated with ADIPOQ genetic variants and lower serum adiponectin levels. Lifestyle Genom. 2020;13(6):172–9. https://doi.org/10.1159/000510021.

    Article  CAS  PubMed  Google Scholar 

  62. Shi R, Ying S, Li Y, Zhu L, Wang X, Jin H. Linking the YTH domain to cancer: the importance of YTH family proteins in epigenetics. Cell Death Dis. 2021;12(4):346. https://doi.org/10.1038/s41419-021-03625-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol. 2019;20(10):608–24. https://doi.org/10.1038/s41580-019-0168-5.

    Article  CAS  PubMed  Google Scholar 

  64. Du H, Zhao Y, He J, Zhang Y, Xi H, Liu M, et al. YTHDF2 destabilizes m(6)A-containing RNA through direct recruitment of the CCR4-NOT deadenylase complex. Nat Commun. 2016;7:12626. https://doi.org/10.1038/ncomms12626.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Zaccara S, Jaffrey SR. A unified model for the function of YTHDF proteins in regulating m6A-modified mRNA. Cell. 2020;181(7):1582–95. https://doi.org/10.1016/j.cell.2020.05.012.

  66. Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, et al. N6-methyladenosine modulates messenger RNA translation efficiency. Cell. 2015;161(6):1388–99. https://doi.org/10.1016/j.cell.2015.05.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Huang H, Weng H, Chen J. m6A modification in coding and non-coding RNAs: roles and therapeutic implications in cancer. Cancer Cell. 2020;37(3):270–88. https://doi.org/10.1016/j.ccell.2020.02.004.

Download references

Acknowledgements

This work was supported by funds from the National Natural Science Foundation of China (Grant No. U21A20249) and China Postdoctoral Science Foundation (2022 M712794). We are deeply grateful to Prof. Lusheng Huang and all the colleagues from State Key Laboratory of Pig Genetic Improvement and Production Technology, Jiangxi Agricultural University, Nanchang for supporting the samples and critical suggestions.

Author information

Author notes

  1. Huanfa Gong and Tao Gong contributed equally to this work.

Authors and Affiliations

  1. Key Laboratory of Molecular Animal Nutrition, Ministry of Education, College of Animal Sciences, Zhejiang University, Hangzhou, 310058, People’s Republic of China

    Huanfa Gong, Tao Gong, Youhua Liu, Yizhen Wang & Xinxia Wang

  2. Key Laboratory of Animal Nutrition and Feed Science in Eastern China, Ministry of Agriculture, College of Animal Sciences, Zhejiang University, Hangzhou, 310058, People’s Republic of China

    Huanfa Gong, Tao Gong, Youhua Liu, Yizhen Wang & Xinxia Wang

Authors
  1. Huanfa Gong
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tao Gong
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Youhua Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yizhen Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xinxia Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YW and HG collected the foundation. YW and XW supervised the project. HG, TG and XW wrote the manuscript. HG analyzed the data and performed visualization. TG performed the experiments and visualization. YL helped in designing the experiments. All the authors have read and agreed the published version of the manuscript.

Corresponding author

Correspondence to Xinxia Wang.

Ethics declarations

Ethics approval and consent to participate

The experimental procedures were in compliance with guidelines of the Committee on Animal Care and Use and Committee on the Ethic of Animal Experiments of Zhejiang University (Hangzhou, China).

Consent for publication

Not applicable.

Competing interests

The author declare no competing interests.

Supplementary Information

Additional file 1: Fig. S1.

RNA expression analysis of MeRIP-seq input data. a and b PCA and heatmap of highly expression genes among 20 and (c and d) 19 samples (excluded L96), respectively. e Volcano plot of RNA differential expression gene (P-adjust < 0.05 and fold change > 1.5), DSPP gene in the dotted box for extremely outlier of the figure (log2foldchange = − 25.3; P-adjust = 1.36E−17)

Additional file 2: Fig. S2.

Establishment of lipogenesis model in vitro. a Oil Red O staining of porcine intramuscular preadipocytes at 0, 2, 4 and 8 d after adipogenic induction, n = 3. b RT-qPCR of ADIPOQ of porcine intramuscular preadipocytes at 0, 2, 4 and 8 d after adipogenic induction, n = 3. c RT-qPCR of PPARγ, CEBPβ and aP2 of porcine intramuscular preadipocytes at 0, 2, 4 and 8 d after adipogenic induction, n = 3. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Additional file 3:

 Table S1. 70 genes in MEdarkturquoise module

Additional file 4:

 Table S2. m6A differential modified genes

Additional file 5:

 Table S3. RNA expression differential genes

Additional file 6:

 Table S4. Methylated and RNA expression co-differential genes

Additional file 7:

 Table S5. KEGG and GO analysis of gene set from Table S1-S4

Additional file 8: Table S6.

m6A modification mutation site of ADIPOQ cDNA in CDS and 3′ UTR

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Gong, H., Gong, T., Liu, Y. et al. Profiling of N6-methyladenosine methylation in porcine longissimus dorsi muscle and unravelling the hub gene ADIPOQ promotes adipogenesis in an m6A-YTHDF1–dependent manner. J Animal Sci Biotechnol 14, 50 (2023). https://doi.org/10.1186/s40104-023-00833-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40104-023-00833-4

Keywords

Journal of Animal Science and Biotechnology

ISSN: 2049-1891