Genome-wide DNA methylation analysis in Chinese Chenghua and Yorkshire pigs

The Chinese Chenghua pig (CHP) is a typical Chinese domestic fatty pig breed with superior meat quality characteristics, while the Yorkshire pig (YP) has the characteristics of fast growth and a high rate of lean meat. Long term natural selection and artificial selection resulted in great phenotypic differences between the two breeds, including growth, development, production performance, meat quality, and coat color. However, genome-wide DNA methylation differences between CHP and YP remain unclear. DNA methylation data were generated for muscle tissues of CHP and YP using reduced representation bisulfite sequencing (RRBS). In this study, a total of 2,416,211 CpG sites were identified. Besides, the genome-wide DNA methylation analysis revealed 722 differentially methylated regions (DMRs) and 466 differentially methylated genes (DMGs) in pairwise CHP vs. YP comparison. Six key genomic regions (Sus scrofa chromosome (SSC)1:253.47–274.23 Mb, SSC6:148.71–169.49 Mb, SSC7:0.25–9.86 Mb, SSC12:43.06–61.49 Mb, SSC14:126.43–140.95 Mb, and SSC18:49.17–54.54 Mb) containing multiple DMRs were identified, and differences of methylation patterns in these regions may be related to phenotypic differences between CHP and YP. Based on the functional analysis of DMGs, 8 DMGs (ADCY1, AGBL4, EXOC2, FUBP3, PAPPA2, PIK3R1, MGMT and MYH8) were considered as important candidate genes associated with muscle development and meat quality traits in pigs. This study explored the difference in meat quality between CHP and YP from the epigenetic point of view, which has important reference significance for the local pork industry and pork food processing.


Background
Epigenetic modifications of the genome can have both short-term and long-term effects on gene expression in different environments [1]. In turn, changes in these expression profiles have implications for multiple traits. DNA methylation was the first discovered epigenetic modification and one of the most thoroughly studied [2]. DNA methylation predominantly occurs at the C-5 position of cytosine in cytosine and guanine dinucleotide (CpG) dinucleotides in mammals [3]. Moreover, DNA methylation is critical for mammalian growth and development [4]. DNA methylation is traditionally regarded as a heritable and stable silence marker, which is essential for X-inactivation [5], silencing of genomic elements such as transposons [6], and genetic imprinting [7]. In addition, variation in DNA methylation involves in a wide range of cellular functions and pathologies [1], and DNA methylation also affects muscle growth and development [3]. Recently, the role of DNA methylation dynamics on skeletal muscle development and disease has been reported [8].
As the main meat source and human medical research model [9], the pig has important research value. Long-term domestication and modern breeding have resulted in both genetic variation and epigenetic modification in different breeds in pigs. Yorkshire pig (YP) is an important commercial pig breed with a high growth rate and lean meat [10]. Chenghua pig (CHP) is a Chinese local breed which is famous for superior meat quality [11]. By contrast, there are significant differences in body composition, muscle, and fat content between Chinese local pigs and commercial pigs [12], especially between CHP and YP [13]. Epigenetic variations, and in particular DNA methylation, might not only influence differences between individuals but also between populations [14]. Hence, DNA methylation might contribute to phenotype variation between pig breeds.
Recently, some studies have explored methylation patterns in different pig breeds and tissues. Choi and colleagues reported the DNA methylome profiles of five different tissues [15]. Zhang and colleagues revealed the epigenetic mechanism of hypoxic adaptation in Tibetan and Yorkshire pigs [16]. Wang and Kadarmideen provided an epigenome-wide DNA methylation map of testis by a genome-wide DNA methylation analysis [17]. However, few studies have investigated the different epigenetic patterns between CHP and YP.
The main objective of this study was to explore the DNA methylation differences between CHP and YP by genome-wide DNA methylation analysis and then identify key genes and candidate epigenetic biomarkers associated with these differences of meat quality traits. We identified the differentially methylation regions (DMRs) and differentially methylation genes (DMGs) of CHP and YP to determine some of the important genomic regions and key genes associated with these phenotypic differences and providing new insights into the epigenetic mechanisms underlying the differences between the two pig breeds.

Summary of RRBS data
Approximately 690.32Gb raw data was generated by RRBS from 48 muscle tissue samples of CHP and YP (approximately 14.38Gb raw data per individual). After quality control, approximately 523.72Gb clean data was obtained (approximately 10.91Gb clean data per individual). Besides, approximately 65% of the reads were mapped to the porcine reference genome (Table 2). Moreover, in all individuals, the density of normalized reads mapped to the proximal and distal regions of the chromosomes was higher than that of reads mapped to other regions. Overall methylated cytosines in the CpG/ CHG/CHH (whereby H can be either A, T, or C) context were 51.39%/0.96%/0.7% in CHP and 52.68%/1.04%/ 0.78% in YP, respectively. Besides, C methylated in an unknown context like CN or CHN (whereby N can be either A, T, G, or C) was observed to be 5.8% in CHP and 5.77 in YP. Figure 1 shows CpG-and non-CpGmethylation sites (CHG, CHH, CN, or CHN) in muscle tissue of CHP and YP.

DMGs identified according to DMRs and functional annotation of DMGs
We annotated 466 DMGs from DMRs identified by comparing CHP vs. YP. Besides, 149 DMGs exhibited higher levels of DNA methylation in CHP than in YP (Table S2), while 317 DMGs exhibited lower levels of DNA methylation in CHP than in YP (Table S3).

Discussion
In this study, we found that there were differences in DNA methylation between CHP and YP. The methylation patterns of CHP may help to explain the epigenetic regulation mechanisms of traits.
Bisulfite sequencing is an ideal and practical technique for studying epigenetic modifications of different species and tissues [18], especially DNA methylation, which can detect the DNA methylation level at each base position of the whole genome. However, genome-wide DNA methylation sequencing with high coverage of the whole genome is required to accurately assess the methylation levels at each base position. Thus, RRBS was used in this study because of its high coverage, small data requirement, low cost, and simple operation. Compared to other studies in pigs [16,17,19], this study used a larger population size. Therefore, RRBS is suitable for detecting DNA methylation differences among breeds in this study.
We observed several interesting GO terms and KEGG pathways associated with muscle metabolism and development. The KEGG pathways of Type II diabetes mellitus (enriched with MAPK10, PRKCE, GCK, MTOR, PIK3R1), cAMP signaling pathway (enriched with ADCY1, ACOX3, ADCY5, ARAP3, PIK3R1, MAPK10, GRIN2B, VIPR2, and VAV2), the GO terms of skeletal muscle acetylcholine-gated channel clustering (enriched with COLQ and DNAJA3) and the cAMP-mediated signaling (enriched with ADCY1, ADCY5, and KSR1) were  Annotated with gene, the percentage of CpG sites or differential methylated regions that overlap with gene promoter, exon, intron, or intergenic; b Annotated within CpG, the percentage of CpG sites or differential methylated regions that overlap with CpG island, CpG shore or other regions Fig. 2 The distribution of differentially methylated regions (DMRs) throughout the whole genome in CHP vs. YP. The purple circle represents the hypomethylated DMRs. The orange triangle represents the hypermethylated DMRs. The color on a chromosome represents the gene density identified in CHP vs. YP. As a major metabolic tissue, metabolic-related pathways and GO terms, including Type II diabetes mellitus, cAMP signaling pathway, skeletal muscle acetylcholine-gated channel clustering, and cAMP-mediated signaling was enriched in this study.
The results indicated that DMGs associated with these metabolic processes show significant differences between CHP and YP. The pork pH has an important relationship with muscle metabolism. In this study, the pH 45min (P = 7.78e-10) and pH 24h (P = 1.24e-4) of CHP were higher than those of YP. Therefore, DMGs involved in muscle metabolism were identified in CHP and YP, which suggested that the difference of pH between the two breeds may be influenced by these pathways and related genes. cAMP signaling pathway is a crucial pathway which regulates pivotal physiologic processes including metabolism, secretion, calcium homeostasis, muscle contraction, cell fate, and gene transcription. In this study, 9 DMGs are enriched in the cAMP signaling pathway. Two of these DMGs, including ADCY1 and PIK3R1, are related to melanoma metastasis. Previous studies have  shown that knockdown of ADCY1 gene leads to decreased intracellular cAMP and subsequently inhibits PKA activity, and phospho-cAMP-responsive element binding protein (CREB) and microphthalmia-associated transcription factor (MITF) levels were significantly downregulated after inactivation of PKA [20]. Furthermore, CREB and MITF have been implicated in melanoma tumor growth and metastasis [21][22][23]. Besides, the ADCY1 gene was identified as a key candidate gene involved in melanoma metastasis [24]. There is an important link between pigmentation and melanoma. This result suggests that ADCY1 gene may affect pigmentation through cAMP. The PI3K protein, encoded by PIK3R1 gene, is a key protein involved in the PI3K/AKT signaling pathway, which is essential for myogenic differentiation [25] and regulates cell survival, growth, differentiation, glucose transport, and utilization [26]. Therefore, the high levels of methylation of ADCY1 and PIK3R1 in CHP may trigger changes in their expression, potentially leading to different meat color traits between CHP and YP. Notably, some other key DMGs, including AGBL4, EXOC2, FUBP3, PAPPA2, MGMT, and MYH8 were found in this study. The AGBL4 gene was regarded as a candidate gene associated with the heterotic quantitative trait in beef cattle [27]. A genome-wide association study (GWAS) suggested that one SNP (rs12210050) in EXOC2 was related to the tanning ability of Europeans [28]. A previous study demonstrated that the FUBP3 gene was associated with the skeletal formation in Duroc population [29]. Furthermore, the FUBP3 gene was identified as a candidate gene associated with the loin eye area in pigs [30]. The PAPPA2 gene encodes pregnancyassociated plasma protein A2 (PAPPA2) which plays an important role in the regulation of IGF-I bioavailability [31]. It is a metalloproteinase that can specifically clew IGFBP-3 and IGFBP-5, thereby releasing IGF-I from its ternary complex, enabling it to bind to IGF-I receptors on the cell surface, initiating growth-promoting activity [32]. Besides, in genome-wide association analysis, PAPP A2 and its related gene, PAPPA, were identified as common genetic variants associated with adult stature in the general population [33]. The MGMT gene is a DNA repair gene responsible for removing alkylation adducts from the O6-position of guanine in DNA. The promoter CpG island hypermethylation associated gene silencing of MGMT is involved in a wide spectrum of human cancers, including glioblastoma [34], gastric [35], colon [36], and ovarian [37]. The MYH8 gene belonged to the myosin heavy chain gene family that share the common features of ATP hydrolysis, actin binding, and potential for kinetic energy transduction [38]. Moreover, the MYH8 myosin is re-expressed during muscle regeneration and is deemed as a specific marker of regenerating fibers in the pathologic skeletal muscle [39,40].

Conclusion
This study performed epigenome-wide DNA methylation analysis using RRBS data generated for muscle tissues of 48 pigs. CHP vs. YP revealed 722 DMRs and 466 DMGs based on these DMRs. Besides, 6 key genomic regions and 8 key DMGs, which might be related to phenotypic differences between CHP and YP, were identified according to the further functional analysis. Our finding may help to further understand the epigenetic mechanisms of phenotype traits and have reference significance for the local pork industry.

Animals and measurements of meat quality
Totals of 48 healthy pigs were used in this study from two pig breeds, including CHP (n = 20) and YP (n = 28). These pigs were maintained in a similar environment to avoid the effects of other confounders. There are 10 males and 10 females in the Chenghua pigs, and there are 20 males and 8 females in the Yorkshire pigs. Each population contains a certain number of males and females. In addition, a large sample size was used to reduce the influence of confounders. Animals were slaughtered at a commercial slaughterhouse when they

Library construction
Briefly, genomic DNA was isolated from flash frozen muscular tissue. Then, the construction of RRBS libraries and paired-end sequencing using Illumina HisSeq analyzer was performed at Novogene technology co., LTD (Beijing, China). Raw sequencing data were processed by an Illumina base-calling pipeline. Genomic DNA was digested with MspI enzyme at 37°C for 16 h. The DNA fragments after enzyme digestion were repaired at the end, and the sequencing adapters with all cytosine methylated were attached. The inserted DNA fragments with the length ranging from 40 to 220 bp were selected for glue cutting. Then, Bisulfite conversion was carried out. After that, the unmethylated C was changed to U (after PCR amplification to T), while the methylated C remained unchanged. Finally, PCR amplification was carried out to obtain the final DNA library. Clean reads were obtained from the raw data after removing reads containing adaptor sequences, unknown, or low-quality bases. The process of quality control was carried out using Trimmomatic software [41]. Quality control was adopted to access the high data quality by (1) removing low-quality reads using a sliding window method (SLIDINGWINDOW: 4:15); (2) removing reads including adaptor sequences (ILLUMINACLIP: adapter.fa: 2:30:7:1: true); (3) removing reads with tail quality lower than 3 or with unknown bases (TRAILING: 3).

Enrichment analysis
Significant GO terms and KEGG pathways were selected after filtering with P < 0.01. R package ggplot2 v3.3.2 was used to visualize the significant GO terms and KEGG pathways for the DMGs associated with DMRs.