Differential DNA methylation analysis across the promoter regions using methylated DNA immunoprecipitation sequencing profiling of porcine loin muscle

Background and Aim: Pork leanness and marbling are among the essential traits of consumer preference. To acquire knowledge about universal epigenetic regulations for improving breed selection, a meta-analysis of methylated DNA immunoprecipitation sequencing (MeDIP-seq) profiling data of mixed loin muscle types was performed in this study. Materials and Methods: MeDIP-seq profiling datasets of longissimus dorsi muscle and psoas major muscles from male and female pigs of Landrace and Tibetan breeds were preprocessed and aligned to the porcine genome. Analysis of differential methylated DNA regions (DMRs) between the breeds was performed by focusing on transcription start sites (TSSs) of known genes (−20,000-3000 bases from TSS). All associated genes were further reviewed for their functions and predicted for transcription factors (TF) possibly associated with their TSSs. Results: When the methylation levels of DMRs in TSS regions of Landrace breed were compared to those of Tibetan breed, 10 DMRs were hypomethylated (Landrace < Tibetan), and 19 DMRs were hypermethylated (Landrace > Tibetan), accordingly (p≤0.001). According to the reviews about gene functions, all associated genes were pieces of evidence for their roles in a variety of muscle and lipid metabolisms. Prediction of the binding TFs revealed the six most abundant binding TFs to such DMRs-associated TSS (p≤0.0001) as follows: ZNF384, Foxd3, IRF1, KLF9, EWSR1-FLI1, HES5, and TFAP2A. Conclusion: Common DMRs-associated TSS between the lean-type and the marbled-type loin muscles were identified in this study. Interestingly, the genes associated with such regions were strongly evidenced for their possible roles on the muscle trait characteristics by which further novel research topics could be focused on them in the future.


Introduction
Pork has been accounted for a large-scale world's meat consumption. Lately, the evolved preference of pork consumers significantly affects its amendment in taste, flavor, juiciness, and tenderness [1,2]. These issues thus become one of the most critical concerns in fields of both veterinary and animal health science as pork leanness and marbling are two essential traits contributing to its quality [1][2][3]. Diverse cascading gene transcription control by epigenetic modification was well-recognized for its involvement in such phenotypic variations among pig breeds [4][5][6][7]. The knowledge in epigenome concerning porcine musculogenesis and lipogenesis has thus significantly contributed to advance in swine breeding for decades [3][4][5][6]8]. For further explanation, epigenetic gene regulation by DNA methylation involves primarily cis-acting DNA elements which are the non-coding DNA binding sites for the gene transcription factors (TFs). Apart from the promoter region, which commonly locates near the transcription start site (TSS), the chromatin loop formation also allows the long-distant cis-acting elements such as enhancers and repressors to cooperate in regulation [9]. Cis-acting elements in mammals, including pigs, usually contain repeated CpG bases, by which their methylation can interfere TF binding. The methylation patterns are inheritable among animal breeds, including swine, and responsible for different transcription patterns among the tissues of interest [4][5][6]9,10]. To study the intricate regulation patterns of genome-wide DNA methylome, large-scale identification of enriched methylated sequences is thus required. Isolation of these methylated DNA fragments could be archived by the antibody against 5-methyl-cytosine, the purification technique recognized as methylated DNA immunoprecipitation (MeDIP) [11]. Coupling of MeDIP and short-read sequencing technologies even allow the high-throughput sequencing of methylated DNA fragments, the technology recognized as MeDIP-sequencing (MeDIP-seq) [12]. Despite its recent introduction, MeDIP-seq has been applied with a few cell types to discover differential DNA methylation patterns among different swine breeds [4,5,7,11,13]. Despite the insight into the epigenetic background of swine breeds implied by the past studies, some limitations in such knowledge's application should also be concerned about.
DNA methylation profiling data of a tissue sample are a composited methylation pattern acquired from common cell types in the particular tissue [14]. Methylation profiling of each tissue type should not be considered as a good representative among one another [14,15]. Fortunately, differential DNA methylation analysis between lean-type and obese-type pig breeds was previously performed using the MeDIP-seq technology. With few replicated muscle tissue samples, a pooled library of one muscle type was constructed entirely for both DNA-seq and mRNA-seq analyses [4]. In the analysis of sequencing data, numbers of discrete biological replicate crucially influence representativeness of the differential analysis outcome [16,17]. By mean of this, the inclusion of numbers of different muscle type biosamples in differential methylation analysis should result in different significant patterns which are universal among various muscle-type tissues.
It should be noted that distances of cis-acting DNA elements from their corresponding TSSs could vary from few to further than 20 Kb [9,18]. While methylated elements are commonly presumed to regulate the transcription of the nearest gene promoter, their distant functions were also evidenced [19,20]. With clusters of different TFs possibly associated with these elements, they can complexly control the transcription of the targeted genes [19,20]. Virtualization of the methylated elements' distances from the TSS and their associated TFs' prediction is thus crucial for productive genome-wide DNA methylation profiling observation. Despite such feasibility, procedures to attain such knowledge in porcine muscle MeDIP-seq study were still limitedly demonstrated.
Gradual increase of porcine meat muscle MeDIP-seq profiling data in the public database has offered us an opportunity to demonstrate the analytical process with them. With available datasets acquired from two types of joint meat muscles, differential DNA methylation profiling analysis between lean and marbled meat pig breeds was performed in this study. Distance determination and TF prediction respectable to the regions were also demonstrated. By including numbers of different muscle tissue datasets, methylated DNA regions associated with the several well-recognized genes between lean and marbled meat muscle were successfully identified.

Ethical approval
All datasets used in this study were available in NCBI SRA public site, and no ethical approval was required.

Sample datasets
The list of porcine meat muscle tissue MeDIP-seq profiling datasets is provided in Table-1. Longissimus dorsi muscle and psoas major muscles from both male and female pigs were included in this study. These muscle tissues were generated from two swine breeds: Landrace as the lean-type breed and Tibetan as the marbled-type meat breed.

Study period and location
The research was conducted at Department of Veterinary Science, Faculty of Veterinary Medicine, Rajamangala University of Technology Tawan-OK, Chon Buri, Thailand, from August 2019 until the end of January 2020.

Analytical workflow
The analytical workflow of this study is shown in Figure-1. The preprocessed sample datasets were aligned to the porcine genome to acquire aligned datasets and counted for the methylated DNA sequence reads (read counts). Differential DNA methylation analysis between muscle datasets of Landrace and Tibetan pigs was performed. Both differential methylation results and read counts in Reads Per Kilobase Millions were visualized by the Circos plot ( Figure-1). Differentially methylated DNA regions (DMRs) between Landrace and Tibetan pigs were subsequently determined for their proximity to TSSs of the associated genes (−20,000-3000 bases from TSS) and displayed by a customized plot (linear genomes displaying in Figure-1). All the associated genes were reviewed for their functions. Finally, TFs associated with these regions would then be predicted ( Figure-1).

Data preprocessing and quality assessment
The sequence read archive (SRA) files of porcine muscle MeDIP-Seq datasets were obtained from the SRA database (https://www.ncbi.nlm.nih.gov/sra) ( Table-1). Data preprocessing was performed as previously described with some modifications [21]. Briefly, each dataset was aligned to the porcine reference genome (Sus scrofa UCSC susScr3) with the Bowtie2 aligner (version 2.3.4.3) [22]. Duplicated sequences of the aligned datasets were then removed [23]. Phred quality score was determined by FastQC [24]. The assessment of MeDIP-Seq quality and DNA methylation score was performed by Bioconductor package MEDIPS (250-bp windows) after extending read to a length of 300 bp [25].

DMRs between landrace and Tibetan pigs
DMRs between Landrace and Tibetan pigs were identified (p≤0.001) by edgeR (exact test for negative binomial distribution) integrated with the MEDIPS package. TSSs of porcine genes were annotated using the porcine genome (Sscrofa 11.1) acquired from the Ensembl database (version 95) (https://asia.ensembl. org/info/data/ftp/index.html). This study only considered the DMRs within −20,000-+3000 bases from TSS. Gene promoter region was considered from −1000 to +500 bases from TSS and was exclusively remarked for the presence of differential methylation within the regions in the karyoplot (Figure-2).

Prediction of binding TFs
RSAT Matrix-Scan (http://rsat.ulb.ac.be/) was used to detect TF binding sites (http://rsat.sb-roscoff. fr/) among the DMRs. Individual sites were predicted using the various position matrix of non-redundant TF binding profiles of vertebrates (JASPAR CORE, non-redundant vertebrates database) and filtered using a threshold on p=10 −6 . For matrix-scan, the first-order Markov chain background model was calculated from the input sequence set.

Virtualization
Count reads of overall methylated DNA regions along with their DMRs between Landrace and Tibetan pigs were illustrated by the Circos plot using the Bioconductor OmicCircos package [26]. To display the DMRs associated with the target genes of interest, the karyoploteR package [27].was utilized to draw the plot for displaying the locations of DMRs and their associated TSSs.

Genome-wide DNA methylation profiling and their differentially methylated DNA regions
After removing the contaminated, low quality, and adapter sequences, each preprocessed dataset had approximately 64 million clean reads (40-42 bps) with at least an 85% alignment rate to the S. scrofa reference genome. The data showed sufficient sequence depth, and coverage analyses revealed that approximately 95% of all CpGs in the S. scrofa genome was covered at least one-fold, whereas nearly 50% of CpGs were covered more than five-fold. Patterns of canonical epigenetic DNA methylation observed among longissimus dorsi muscle and psoas major muscle tissues acquired from both Landrace and Tibetan breeds were similar, by which strong positive correlation across the breed was demonstrated by correlation analysis (r=0.75-0.95, p<10 −16 ). The comparison of the DNA methylation profiles of the Landrace and Tibetan breeds revealed 74,009 DMRs (p≤0.001) Available at www.veterinaryworld.org/Vol.13/June-2020/14.pdf      d after merging significant neighboring windows. When compared, methylation levels of Landrace breed in these DMRs to those of Tibetan breed, 28,028 hypomethylated regions (Landrace < Tibetan) and 45,981 hypermethylated regions (Landrace > Tibetan) were found scattering among somatic and X chromosomes (Figure-3).

Discussion
DNA methylation patterns mutually distributed among the tissue types usually implied universal metabolic features shared in common among them [14]. Based on such an approach, this study included considerable MeDIPS datasets generated from two common loin meat muscle tissues, longissimus dorsi muscle and psoas major, with the Available at www.veterinaryworld.org/Vol.13/June-2020/14.pdf differential methylation analysis to discover the universal differences in DNA methylation pattern of loin muscle to understand the differences between lean and marbled meat pig breeds. Since only rigid cellular metabolisms were controlled in universal patterns among tissue types, limited numbers of genes with differentially methylated TSS regions were expected as successfully demonstrated in this study.
Several novel genes with TSS-associated DMRs were also exclusively identified in this study. Apart from constructing sample libraries from different meat tissues, other factors, such as muscle types, dataset numbers, and pig breeds, were also considered as contributors to such achievement (Figure-2). For the available virtualization of these regions' localization following their associated genes, the karyoplot was successfully applied [27]. Among all the presented TSS-associated DMRs, some DMRs in the promoter regions of three coding genes, ESSSCG00000038672 (hCH-CHD2), ESSSCG00000027257 (PSMB1), and ESSSCG00000031489 (DUSP22) were also

Chr
GeneID Common name Function References chr1 ENSSSCG00000027257 PSMB1 Upregulated according to obesity, PSMB1 was transcriptional activator of RBP4 -a gene associated with insulin resistance and transcription activation of adipocyte. [28,29] chr2 ENSSSCG00000014560 COX8H COX8H was related to cytochrome C oxidase -the crucial enzyme with potential role in oxidative and fatty acid metabolisms of muscle fibers. [30,31] ENSSSCG00000014561 NLRP6 NLRP6 was the sensor component of NLPP3 inflammasome which could enhance muscular dystrophy and was strongly associated with obesity. [7,[32][33][34] ENSSSCG00000025023 PGGHG PGGHG was a catalyzing enzymes of D-glucose, and thus could affect both muscle and adipose cell's metabolisms. [35,36] chr3 ENSSSCG00000007542 PRKAR1B PRKAR1B encoded a regulatory subunit of cyclic AMPdependent protein kinase A in the signaling pathway of the second messenger cAMP -contributing to myofiber size and negatively associating with intramuscular fat. [32,37] chr6 ENSSSCG00000026662 hPRDM7 PRDM7 could function as a histone methyltransferase, and thus was likely affect DNA methylation affecting muscle tissue development and function. [38] ENSSSCG00000029761 hURAD URAD enzyme catalyzed anti-oxidant process of uric acid to form allantoin -the marker of oxidative stress in muscle. [39,40] chr7 ENSSSCG00000031489 DUSP22 DUSP22 protein activate the JNK signaling pathwayregulating the phosphorylation state of several kinases in skeletal muscle and requiring for obesity.
[ [45][46][47] chr10 ENSSSCG00000010801 CDC73 CDC73 protein was a subunit of PAF protein complex associated with the RNA polymerase II subunit and a histone methyltransferase complex, and thus played role in chromatin modifications and gene transcription of muscle tissue. [48] ENSSSCG00000031991 GLRX2 GLRX2 protein implicated MITOCHONDRIAL REDOX REGULATION preventing oxidative damage and was evidenced to control oxidative phosphorylation in cardiac muscles. [49,50] chr11 ENSSSCG00000022638 ATP12A ATP12A was responsible for potassium absorption for muscle contraction -by which the lack of potassium could result in muscle destruction by decreased blood flow. [51,52] chr12 ENSSSCG00000017137 METRNL METRNL product was associated with adipocyte browning and glucose tolerance. [53] ENSSSCG00000027229 Several miRNAs in skeletal muscle could regulate myogenesis, hypertrophy, atrophy, and regeneration.
[54] chr13 ENSSSCG00000011178 CPNE4 CPNE4 product was calcium-dependent phospholipidbinding protein, and thus can involve in several calciummediated processes in muscle tissue.
[55] chr14 ENSSSCG00000038672 hCHCHD2 CHCHD2 was a transcription factor of cytochrome C oxidase required for optimal mitochondrial function of muscle cell. [56] chr15 ENSSSCG00000022011 NMI N-myc can functionally replace c-myc in murine development, cellular growth, and differentiation. [57] chr16 ENSSSCG00000032942 DAP DAP was necessary for complete myotube development during muscle fiber formation.
[58] chrX ENSSSCG00000038044 GYG2 Glycogenin was important for muscle glycogenesis -by which glycogen over-accumulation could result in fatigue reduced function of muscle. [59] Available at www.veterinaryworld.org/Vol.13/June-2020/14.pdf noticeable ( Figure-2a-c). Due to the robust influence of promoter on epigenetic gene regulation [9,19]. these differentially methylated genes should be considered as exclusive candidates for further transcription determination. Of note, the current study only focused on TSS-associated methylated regions. Further analysis of other potential epigenetic gene regulation regions, including gene exons [60] and introns [61,62], was thus encouraged in separated studies. In this study, several novel genes in loin muscles manifested their differences in the TSS region methylation levels in Landrace comparing to Tibetan pigs. Landrace pig is recognized for its lean large muscle bundle, while Tibetan pig muscle is more compact with outstanding marbling levels [63]. The different characteristics of the two breeds thus greatly denoted the diverge metabolisms of their muscle tissues. Supporting this, almost all differentially methylated genes also manifested their evidenced functions relating to musculogenesis and lipid metabolisms in this study (Table-2). Interestingly, two differentially methylated miRNA genes, ENSSSCG00000020190 and ENSSSCG00000020232, were also identified ( Figure-2d), contributing new targets for a miRNA study in the skeletal muscles [54].
The regulatory function of the acquired TSSassociated DMRs required further biological evaluation despite the knowledge of their methylation levels. Since DNA methylation could either inhibit or recruit particular TFs in various circumstances [9,20,64]. the prediction of candidate TFs binding regions was thus crucial to study their epigenetic gene regulation. Due to the limited demonstration of such process in porcine skeletal muscle MeDIP-seq study, we introduced the RSAT Matrix-Scan as a convenient web-based tool for the prediction of TF binding. Of note, this study intentionally showed only TFs with the highest prediction significant score for demonstration. Since each TSS region could be occupied by various TFs [65]. other significantly predicted TFs (p≤10 −6 ) should also be considered for any further intensive study.

Conclusion
In summary, the determination of TSS regions that were differentially methylated between lean-type and marbled-type pig breeds in loin meat muscle tissues was demonstrated as a novel approach in this study. The procedure rendered the limited numbers of differentially methylated genes, but with strong concatenation with the different muscle characteristics presented between these two breeds. We also introduced a karyoplot for the available virtualization of TSS-associated DMRs along with their candidate TFs binding to them. As a result, the knowledge acquired from this study established several new research topics in lean and marbled meats to evaluate the epigenetic gene regulation. It is worth noting that this study aimed to demonstrate a beneficial strategy for porcine muscle MeDIP-Seq data analysis, which was also applicable to other biological systems.

Authors' Contributions
KC, TS and RK planned the study design, collected the datasets, and analyzed the data. EO and CN refined the study design and the objective. DP carried out technical coding correction and hardware maintenance. KC and RK drafted and reviewed the manuscript. All authors read and approved the final manuscript.