Identifying key DNA methylation sites and their cis-methylation quantitative loci for intramuscular fatty acid traits using genome and methylome data in Yorkshire pigs

Abstract

Single-nucleotide polymorphisms (SNPs) that are located near CpG sites have the potential to influence DNA methylation levels, leading to distinct patterns of DNA methylation. These methylation quantitative trait loci (meQTLs) can influence methylation levels across the entire genome, thereby impacting phenotypic traits. In this study, we utilised reduced representation bisulphite sequencing (RRBS) to analyse DNA methylation in muscle tissue samples collected from 140 Yorkshire pigs. Additionally, we assessed SNPs using the Porcine 50K chip, examined 17 fatty acid traits, and conducted epigenome-wide association studies (EWAS) and genome-wide association studies (GWAS). Our investigation revealed numerous associations between DNA methylation and fatty acid traits, including C16:0, C18:0, C20:3n-3, and C22:6n-3. Notably, a substantial portion of these associations was specifically identified by EWAS and not GWAS. Subsequently, we performed GWAS to identify meQTLs using methylation levels of significant CpGs from the EWAS as traits. The results revealed the predominance of trans-meQTLs; however, we also discovered a cis-meQTL on SSC12. Finally, we reported a crucial CpG site (SSC12:3735523) and its cis-meQTL on SSC12, which are located in close proximity to the FASN gene. Our findings contribute to the identification of CpG sites and their cis-meQTLs associated with fatty acid traits, utilising comprehensive genome and methylome data. Moreover, our study provides novel insights into the molecular mechanisms underlying fatty acid traits in pigs.

Keywords:

• Pig
• DNA methylation
• EWAS
• GWAS
• meQTL

Introduction

The intramuscular fatty acid composition, serving as the fundamental energy source and a crucial regulator of cellular metabolism in animals, holds significant importance in terms of its nutritional properties and impact on human health (Wood et al. 2004). Fatty acids encompass various forms, including saturated fatty acids (SFAs), trans-fatty acids (TFAs), monounsaturated fatty acids (MUFAs), and polyunsaturated fatty acids (PUFAs) (Kremmyda et al. 2011; Tvrzicka et al. 2011). Additionally, the fatty acid composition substantially influences meat quality traits. Elevated levels of PUFAs have been linked to diminished oxidative stability and flavour (Wood et al. 2008), while augmenting the content of MUFAs has shown potential in enhancing organoleptic and technological qualities, as well as nutritional properties (Cameron et al. 2000). Therefore, acquiring a comprehensive understanding of the underlying mechanisms that govern pig fatty acid composition is crucial for enhancing the quality of pork and optimising production efficiency.

In the past few decades, extensive genome-wide association studies (GWAS) conducted in diverse pig populations have successfully identified numerous genetic loci associated with fatty acid traits (Nii et al. 2006; Guo et al. 2009; Muñoz et al. 2013; Ros-Freixedes et al. 2016; Zhang et al. 2016). Additionally, candidate genes involved in lipid metabolism, such as ELOVL fatty acid elongase 6 and 7 (ELOVL6 and ELOVL7) (Zhang et al. 2016), fatty acid synthase (FASN) (Sato et al. 2017), fatty acid desaturase 1 and 2 (FADS1 and FADS2) (He et al. 2018), and stearoyl-CoA desaturase (SCD) (Viterbo et al. 2018) have been discovered for fatty acid traits in adipose and muscle tissue. Similarly, (Orozco et al. 2014) highlighted the inter-individual variability exhibited by epigenetic modifications, much like genetic variations. Methylation quantitative trait locus (meQTL) studies in human whole blood have identified approximately 4.7 million cis-meQTLs and 630 thousand trans-meQTLs targeting over 120 thousand CpG sites (Huan et al. 2019). These findings are supported by (Verma 2012; Xu et al. 2018; Gomez-Alonso et al. 2021). (Schulz et al. 2017) analysed 110 human hippocampal tissues, and they discovered that DNA methylation in hippocampal tissues is linked to various complex traits and diseases, resulting in the identification of cis-meQTLs at 14,118 CpG loci. They further integrated these findings with gene expression data to construct a comprehensive map of meQTLs and expression quantitative trait loci (eQTLs) in the human hippocampus. In a separate study, Lin et al.(Lin et al. 2018) explored meQTLs in 344 human post-mortem brain tissue samples, constructing a genetic regulatory map of DNA methylation in the human brain genome. Their study revealed that approximately 55% of CpG loci across the genome exhibited meQTL associations, indicating the widespread presence of meQTLs throughout the genome. Furthermore, (Gong et al. 2019) in their study on various cancer tissue samples, observed a higher proportion of trans-meQTLs. They identified trans-meQTLs in more than 10% of the total CpG loci analysed in relation to cancer. Collectively, these studies provide further evidence supporting the importance of identifying meQTLs in enhancing our understanding of the impact of SNPs on CpG loci methylation.

In light of the hypothesis that genetic variation can influence DNA methylation and subsequently impact fatty acid traits in pigs, our study aimed to elucidate the molecular mechanisms underlying the genetic regulation of porcine fatty acid traits. We employed an integrated approach that combined genomic and methylomic data to achieve this. By doing so, we aimed to identify crucial DNA methylation sites and their corresponding cis-meQTLs associated with fatty acid traits in pigs.

In our investigation, we conducted epigenome-wide association studies (EWAS) to elucidate the intricate relationship between DNA methylation levels and 17 fatty acid traits in a carefully selected cohort of 140 Yorkshire pigs. The primary objective of our study was to determine the extent to which DNA methylation contributes to the observed variation in fatty acid traits. To accomplish this, we integrated SNP data with reduced representation bisulphite sequencing (RRBS) data, enabling us to reveal a substantial number of associations between DNA methylation patterns and fatty acid traits. Notably, a significant proportion of these associations remained undetected when employing conventional genome-wide association studies (GWAS) methods. Furthermore, to gain further insights into the genetic underpinnings of the observed associations, we conducted GWAS using the significant CpG sites identified from the EWAS as phenotypes and the SNPs as predictors. Through this comprehensive analysis, we successfully identified a cis-meQTL associated with a key CpG site situated at position SSC12:3735523 on the pig genome. Remarkably, the FASN gene, which plays a crucial role in fatty acid metabolism, was found to be in close proximity to the identified key CpG site and its corresponding cis-meQTL. These compelling findings provide compelling evidence of the involvement of DNA methylation in modulating fatty acid traits in pigs and offer valuable insights into the molecular mechanisms underlying the role of the FASN gene in fatty acid metabolism.

Material and methods

Ethics declarations

All experimental procedures were conducted in strict compliance with the ethical guidelines and regulations set forth by the Institutional Review Board (IRB14044) and the Institutional Animal Care and Use Committee of the Sichuan Agricultural University, with complete adherence to permit number DKY-B20140302. These measures ensured the ethical treatment and welfare of the animals involved in the study.

Animals and DNA sample

A total of 140 Yorkshire pigs consisting of 51 males and 89 females, were utilised to extract genomic DNA from muscular tissue. Each pig was considered to be an individual experimental unit. The pigs were raised under the same recommended environment at the conservation farm of Mingxing Agricultural Science and Technology Development Co., Ltd. The pigs with an initial weight of about 30 kg were group-housed in cement-floor pens with 20 pigs per pen. The pigs were fed twice daily with the same diet, and water was provided ad libitum using nipple drinkers. The diet composition was formulated using corn and soybean meal, ensuring that it contained adequate concentrations of crude protein, trace minerals, and vitamins that met or exceeded the recommendations set by the National Research Council (NRC, 1998) for different stages of growth (Luo et al. 2021). All individuals were raised until they reached an average weight of 111.71 kg ($\pm 12.97$ kg) before being transported to the slaughterhouse. A 24-h fasting period was observed prior to slaughter determination. Once the carcase composition traits were determined, samples were taken from the region between the fourth and fifth last rib for fatty acid composition analysis. Samples from the muscle tissue were snap-frozen in liquid nitrogen and stored at −80 °C until analysis. The fatty acid composition was analysed using a GC-¹⁴C gas chromatograph (Shimadzu, Kyoto, Japan) following a previously published method (Qin et al. 2007). The fatty acid components in the test sample were determined by comparing the peak times of each component in the fatty acid standard solution using gas chromatography. The percentage content of each fatty acid component was calculated using the peak area normalisation method.

Library construction

Genomic DNA was extracted from flash-frozen muscular tissue samples. Subsequently, RRBS libraries were constructed, and paired-end sequencing was conducted using the Illumina HiSeq analyser at Novogene Technology Co., Ltd. (Beijing, China). The raw sequencing data underwent processing using the Illumina base-calling pipeline. The genomic DNA was first digested with the MspI enzyme at 37 °C for 16 h. The resulting DNA fragments were repaired at the ends, and sequencing adapters with all cytosines methylated were ligated. DNA fragments with a length ranging from 40 to 220 bp were selected for size selection. Following this, bisulphite conversion was performed. During bisulphite conversion, unmethylated cytosines were converted to uracils (which were subsequently converted to thymine during PCR amplification), while methylated cytosines remained unchanged. Finally, PCR amplification was conducted to obtain the final DNA library. The Trimmomatic software (Bolger et al. 2014) was used to obtain clean reads from the raw data. The following quality control criteria were applied: (1) removal of low-quality reads using a sliding window method (SLIDINGWINDOW:4:15); (2) removal of reads containing adapter sequences (ILLUMINACLIP:adapter.fa:2:30:7:1:true); and (3) removal of reads with a trailing quality score lower than 3 or with unknown bases (TRAILING:3).

SNP data

The genotyping of pigs was performed using the Neogen GeneSeek Porcine50k chip, which is based on the Illumina platform (Illumina Inc., San Diego, USA). Genotypes were subjected to a series of filtering criteria. Loci were filtered to ensure they met the following criteria: Hardy-Weinberg equilibrium $P>1.0\times {10}^{-6},$ minor allele frequency (MAF) > 0.05, and call rate > 0.9. After applying these quality control measures, the final dataset used for subsequent analyses comprised a total of 36,163 loci.

Data analysis

The clean reads obtained from sequencing were aligned to the pig reference genome (Sscrofa11.1) using the Bismark software (Krueger and Andrews 2011) after removing adaptor sequences. Subsequently, the methylation levels of CpG sites were measured using the bismark_methylation_extractor program. To ensure high-quality data, we filtered the methylation level data by retaining only CpG cytosines that were present in all samples and had coverage in at least 30 samples. Furthermore, we applied an additional filter to include CpG sites with a minimum coverage of 10 reads within each sample.

The distribution of CpG sites was analysed using the R package Rldeogram v0.2.2 (Hao et al. 2020). The R package genomation v1.20 (Akalin et al. 2015) was used to annotate all the analysed CpGs. The annotations were based on the porcine RefSeq and CpG island databases (Sscrofa11.1/susScr11), obtained from the UCSC website (http://genome.ucsc.edu/cgi-bin/hgTables). In our analysis, the gene body region was defined as the region from the transcription start site (TSS) to the transcription termination site, while the promoter region was defined as the 2-kb region upstream of the TSS.

Epigenome-wide association study

The R package CpGassoc v2.60 was used to assess the association between the methylation level of CpG sites and 17 fatty acid traits. To perform this analysis, we employed a linear mixed-model approach, considering the methylation level of CpG as the outcome variable. The model was adjusted for covariates, including sex, farm, and days of age. The specific model applied was as follows: (1) $$y=\mu +x\beta +u+e$$(1) where $y$ is the vector of the phenotypes, $\mu$ is the mean, $x$ is the vector of CpG, $\beta$ is the CpG effect, $u$ represents random effects, $e$ is the vector of random residuals.

We corrected for multiple comparisons using a Bonferroni correction for 0.05/4863340 test setting a significant threshold P value to $1.0821 \times {10}^{-8}.$ Subsequently, we constructed Manhattan plots to present the results of epigenome-wide association analysis. The Manhattan plots and quantile-quantile plots (Q-Q plots) were performed using the R package qqman v0.1.4 (Turner 2014).

Candidate genes and annotation

To identify candidate genes associated with the significant CpG sites, we utilised the Ensembl Biomart database (http://www.biomart.org). For intragenic associations, the candidate gene corresponds directly to the gene itself. For intergenic associations, we considered the two nearest flanking genes based on their genomic distance from the CpG site. The distances between the site and each flanking gene were recorded for intergenic associations.

To further investigate the functional implications of these candidate genes, we performed enrichment analysis using the Database for Annotation, Visualisation, and Integrated Discovery (DAVID) v6.8 (http://david.ncifcrf.gov/). The gene lists obtained from the previous step were submitted to DAVID for analysis. Significant Gene Ontology (GO) terms and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways were identified, filtering the results with a significance threshold of p < 0.01.

Genome-wide association study for fatty acid traits

GWAS analyses were conducted independently for each of the 17 fatty acid traits using the GEMMA software (DePristo et al. 2011). The association analysis was performed using a mixed linear model, which accounted for the complex structure of the data. In this model, the phenotype traits were considered as the dependent variable, and adjustments were made for potential confounding factors such as sex, farm, and days of age. The specific model used for the analysis can be described as follows: (2) $$\mathrm{y}=\mathrm{Xm}+\mathrm{Wa}+\mathrm{u}+\mathrm{e}$$(2) where $\mathrm{y}$ is the vector of the phenotypes; $\mathrm{m}$ is the SNP effects; $\mathrm{a}$ is the vector of residual polygene effects; $\mathrm{u}$ represent the random effects; $\mathrm{X}$ and $\mathrm{W}$ are the incidence matrixes for $\mathrm{m}$ and $\mathrm{a},$ respectively; $\mathrm{e}$ is the vector of random residuals. SNP associations that achieved the Bonferroni threshold of $P<1.38\times {10}^{-6}$ were considered as significant associations.

Overlap of EWAS and GWAS

A proximity criterion was applied to determine the overlap between associations obtained from the EWAS and GWAS analyses. Specifically, associations were considered to be overlapping if they were located within a genomic distance of 2 megabases (2 Mb) from each other (Battram et al. 2022). This approach allowed for the identification of potential shared genetic variants or genomic regions that contribute to both DNA methylation patterns and fatty acid traits.

Identification of methylation quantitative trait loci

The methylation levels of CpG loci that remained significant in the previous conditional association analysis were used as dependent variables with SNP as the independent variable and subjected to association analysis to identify methylated quantitative trait loci (meQTL). A mixed linear model was used in the association analysis as follows: $$y=Xm+Wa+e$$ where y is the CpG methylation level; m is the SNP effect; a is the residual polygenic effect; e is the residual; and X and W are the association matrix of m and a, respectively. SNP associations achieve the Bonferroni threshold of $P<4.56\times {10}^{-9}$ were considered as significant associations.

Identification of meQTL by methylation GWAS

GWAS was conducted to identify methylation quantitative trait loci (meQTLs) using the significant CpG sites identified from the EWAS as traits and single nucleotide polymorphisms (SNPs) as predictors. The GEMMA software (DePristo et al. 2011) was utilised for the association analysis, employing a mixed linear model. It is important to note that there is a distinction between the GWAS model described previously, which focused on fatty acid traits with the vector of phenotypes denoted as y, and the methylation GWAS linear mixed model where the vector of methylation levels is represented as y. This difference reflects the specific nature of the analysis, which aims to elucidate the genetic factors influencing DNA methylation patterns.

Results

Phenotype statistics

This study conducted a comprehensive analysis of fatty acid composition on muscle samples obtained from 140 pigs. Table 1 presents the results, revealing the abundance of 17 different fatty acid traits in the samples. Among them, oleic acid (C18:1n-9) was found to be the most predominant fatty acid, constituting over 30% of the total fatty acid content. Following oleic acid, palmitic acid (C16:0) and linoleic acid (C18:2n-6) were observed as the next most abundant fatty acids. Notably, these 17 fatty acid traits collectively account for more than 95% of the total fatty acid content in the muscle samples.

Table 1. Summary of fatty acid Profile of the Yorkshire pigs (% of total fatty acids).

Fat acid traits	Yorkshire pigs (n = 140)
Decanoic acid (C10:0)	0.19$\pm$0.18
Dodecanoic acid (C12:0)	0.10$\pm$0.04
Myristic acid (C14:0)	1.13$\pm$0.19
Palmitic acid (C16:0)	24.10$\pm$1.63
Stearic acid (C18:0)	10.84$\pm$2.27
Arachidic acid (C20:0)	0.14$\pm$0.03
Palmitoleic acid (C16:1)	2.72$\pm$0.53
Oleic acid (C18:1n-9)	33.62$\pm$4.07
cis-11-eicosenoic acid (C20:1n-9)	0.56$\pm$0.13
Linoleic acid (C18:2n-6)	19.03$\pm$4.04
α-linolenic acid (C18:3n-3)	0.70$\pm$0.17
γ-linolenic acid (C18:3n-6)	0.16$\pm$0.05
cis-11,14-eicosadienoic acid (C20:2)	0.64$\pm$0.12
cis-11,14,17-eicosatrienoic acid (C20:3n-3)	0.11$\pm$0.03
cis-8,11,14-eicosatrienoic acid (C20:3n-6)	0.44$\pm$0.13
Arachidonic acid (C20:4n-6)	0.14$\pm$0.16
Docosahexaenoic acid (C22:6n-3)	0.33$\pm$0.18
Saturated Fatty acids (SFAs)^a	39.75$\pm$2.43
Monounsaturated fatty acids (MUFAs)^b	37.87$\pm$4.34
Polyunsaturated fatty acids (PUFAs)^c	21.59$\pm$4.26

Results are presented as the means $\pm$ SD.

Saturated Fatty Acids (SFAs)a = C10:0 + C12:0 + C14:0 +  C16:0 + C18:0 + C20:0 + C21:0+ C22:0+ C23:0;.

Monounsaturated Fatty Acids (MUFAs)b = C14:1 + C15:1 + C16:1 + C17:1 + C18:1n-9 + C20:1n-9 + C22:1n-9 + C24:1n-9;.

Polyunsaturated Fatty Acids (PUFAs)c = C18:2n-6 + C20:2 + C18:3n-3 + C18:3n-6 + C20:3n-3 + C20:3n-6 + C20:4n-6 + C20:5n-3 + C22:2 + C22:6n-3.

Summary of RRBS data

RRBS libraries were constructed from muscular tissue samples obtained from 140 Yorkshire pigs (Table 2), and the sequencing was performed using the Illumina HiSeq platform. On average, each sample yielded approximately 14.02 ± 2.45 Gb of raw data. After quality control, an average of 11.13 ± 2.14 Gb of clean data per sample was obtained. The bisulphite conversion efficiency of the Yorkshire pigs exceeded 99%, indicating a robust conversion of unmethylated cytosines to uracils. Furthermore, more than 60% of the reads from the Yorkshire pigs were successfully mapped to the porcine reference genome.

Table 2. Mapping results of reduced representation bisulphite sequencing (RRBS) data in the Yorkshire pigs.

Items	Yorkshire pigs (n = 140)
Mean raw data (Gb)	14.02$\pm$2.45
Mean clean data (Gb)	11.13$\pm$2.14
Mean bisulphite conversion rate (%)	99.57$\pm$0.07
Mean mapping rate (%)	65.40$\pm$5.35

For further analysis, the methylation data were filtered to include only CpG sites with a minimum coverage of 10× and present in at least 30 samples. This filtering resulted in a total of 4,863,340 CpG sites for subsequent analysis (Table 3). The distribution of these CpG sites across different genomic regions was as follows: 3.73% in the promoter region, 4.67% in exons, 38.41% in introns, and 53.19% in intergenic regions. Additionally, CpG islands (CGIs) are generically equipped to influence local chromatin structure and simplify the regulation of gene activity. Furthermore, there is a strong correlation between CGIs and transcription initiation. The percentages of these CpG sites annotated within CpG islands, CpG island shores, and other regions were 3.20%, 11.54%, and 85.26%, respectively.

Table 3. Annotation of Global CpG sites and significant CpG sites identified for EWAS.

Genetic features	Number	Annotated with gene				Annotated within CpG
		Promoter	Exon	Intron	Intergenic	CpG island	CpG shore	Other regions
Global CpG sites	4863340	3.73%	4.67%	38.41%	53.19%	3.20%	11.54%	85.26%
Significant CpG sites	2390	5.02%	5.31%	42.93%	46.74%	2.59%	10.33%	87.08%

Annotated with gene = the percentage of CpG sites that overlap with gene promoter, exon, intron, or intergenic;.

Annotated within CpG = the percentage of CpG sites that overlap with CpG island, CpG shore or other regions.

EWAS identifies methylation signatures significantly associated with fat acid traits

We conducted an Epigenome-Wide Association Study (EWAS) to investigate the relationship between CpG methylation levels and 17 fatty acid traits. The results of the EWAS analysis are presented in Figure 1.

Figure 1. Manhattan plots of fatty acid traits in Epigenome-Wide Association Study (EWAS). (A) C10:0; (B) C12:0; (C) C14:0; (D) C16:0; (E) C16:1; (F) C18:0; (G) C18:1n-9; (H) C18:2n-6; (I) C18:3n-3; (J) C18:3n-6; (K) C20:0; (L) C20:1n-9; (M) C20:2; (N) C20:3n-3; (O) C20:3n-6; (P) C20:4n-6; (Q) C22:6n-3.

Note: The x-axis represents the chromosome and the y-axis represents the significant of all CpG sites, represented by -log10(P-value). The red line represents the genome level significance threshold based on 0.05/N, where N is the number of CpG sites.

A total of 2,390 significant associations were identified between CpG sites and fatty acid traits (Table 4). Among these significant CpG sites, 5.02% were annotated within promoter regions, 5.31% were annotated within exon regions, 42.93% were annotated within intron regions, and 46.74% were annotated within intergenic regions. Regarding the distribution of these significant CpG sites, 2.59% were annotated within CpG islands, 10.33% were annotated within CpG island shores, and 87.08% were annotated within other regions (Table 3).

Table 4. Summary of Epigenome-Wide Association Study (EWAS) results and Genome-Wide Association Study (GWAS) results for 17 fatty acid traits.

Fat acid traits	EWAS associations	GWAS associations
Decanoic acid (C10:0)	140	0
Dodecanoic acid (C12:0)	13	2
Myristic acid (C14:0)	0	0
Palmitic acid (C16:0)	330	0
Stearic acid (C18:0)	282	0
Arachidic acid (C20:0)	1	0
Palmitoleic acid (C16:1)	0	0
Oleic acid (C18:1n-9)	0	0
cis-11-eicosenoic acid (C20:1n-9)	25	0
Linoleic acid (C18:2n-6)	61	0
α-linolenic acid (C18:3n-3)	0	0
γ-linolenic acid (C18:3n-6)	0	0
cis-11,14-eicosadienoic acid (C20:2)	0	0
cis-11,14,17-eicosatrienoic acid (C20:3n-3)	271	4
cis-8,11,14-eicosatrienoic acid (C20:3n-6)	1	0
Arachidonic acid (C20:4n-6)	1	1
Docosahexaenoic acid (C22:6n-3)	1265	9
Total	2390	16

Based on the 2,390 significant CpG sites, we identified a total of 763 protein-coding genes and lncRNAs (Table S1). These genes and lncRNAs were then subjected to enrichment analysis, resulting in the identification of 18 Gene Ontology (GO) terms and 7 Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways that showed enrichment. Notable enriched terms and pathways included Insulin resistance, actin filament bundle assembly, and F-actin capping protein complex (Figure 2). It is indicated that these genes and lncRNAs may affect meat quality traits by regulating biological pathways related to insulin and actin composition.

Figure 2. Bubble diagram of gene Ontology (GO) terms and Kyoto encyclopaedia of genes and genomes (KEGG) pathway terms for significant associations in Epigenome-Wide Association Study (EWAS).

The X-axis represents the P value of genes enriched in the corresponding GO and KEGG pathway terms. The Y-axis represents the GO and KEGG pathway terms to which the genes are enriched. The bubble shape represents the classification of GO and KEGG pathway terms. The colour of the bubble represents the log transformation of the P value.

Additionally, we used the inflation factor (λ), which represents the regression of observed values on expected values, to assess population stratification. An average inflation factor of 0.95 indicated that there was no inflation in the EWAS results, suggesting that the influence of population stratification was effectively controlled for in the analysis.

Genome-wide association analyses for fatty acid traits

After quality control, a total of 36,163 SNPs were retained for GWAS analysis of the 17 fatty acid traits (Table 4). Among these SNPs, we identified 16 significant associations ($P<1.38\times {10}^{-6}$)with the fatty acid traits. Specifically, there were 2 significant SNPs associated with C12:0, 4 significant SNPs associated with C20:3n-3, 1 significant SNP associated with C20:4n-6, and 9 significant SNPs associated with C22:6n-3. For the C12:0 trait, the two significant SNPs were located on SSC2 at positions 28,709,127 bp and 28,740,723 bp. Regarding the C20:3n-3 trait, the four significant SNPs were located on SSC12 at positions 25,708,148 bp and 43,997,932 bp, as well as on SSC16 at positions 75,205,458 bp and 75,391,932 bp. The significant SNP associated with C20:4n-6 was located on SSC14 at position 12,857,576 bp. Lastly, for the C22:6n-3 trait, the nine significant SNPs were located on SSC4 at positions 106,382,235 bp, 106,485,236 bp, 109,093,503 bp, 109,885,605 bp, and 110,552,282 bp, as well as on SSC12 at positions 60,260,816 bp and 60,372,315 bp, and on SSC13 at positions 94,960,811 bp and 95,481,059 bp.

Identification of the overlap between EWAS and GWAS

In our study, we identified a significant number of associations between fatty acid traits and CpG methylation using EWAS. Importantly, some of these associations were not identified using traditional GWAS methods. To assess the overlap between EWAS and GWAS, we considered associations to be overlapping if the associated loci were within a 2 Mb distance from each other. Among the fatty acid traits studied, we identified a total of 2,390 significant associations through EWAS and 16 significant associations through GWAS. Out of these, 9 associations were found to be overlapping (Figure 3), where the associated loci were within 2 Mb of each other (as indicated in Table S2). This suggests that there are shared genetic and epigenetic factors contributing to the observed associations between CpG methylation and fatty acid traits in our study population.

Figure 3. venn plot showing the overlap between Epigenome-Wide Association Study (EWAS) significant associations and Genome-Wide Association Study (GWAS) significant associations.

The red circle represents the EWAS significant associations, the green circle represents the GWAS significant associations.

Conditional association studies

We conducted conditional association studies to investigate whether the associations identified through EWAS and GWAS were attributable to the same underlying signal or represented co-localizing but independent signals. To achieve this, we performed EWAS for the overlapping associations while considering the SNP genotype as a covariate. In this analytical approach, the disappearance of an EWAS association when incorporating the SNP genotype as a covariate indicates that the overlapping association between EWAS and GWAS resulted from the same underlying signal. Conversely, if the EWAS association remains statistically significant even after adjusting for the SNP genotype, it suggests that the overlapping association is independent but co-localizing.

Among the nine overlapping associations observed (as presented in Table S3), none reached statistical significance at the Bonferroni-corrected threshold ($P<5.56\times {10}^{-3}$). These findings imply that the same signal likely drove the associations we identified using both EWAS and GWAS at the associated locus.

Out of the 9 overlapping associations (Table S3), none were significant at the Bonferroni threshold ($P<5.56\times {10}^{-3}$). The result suggests that the associations we found using both EWAS and GWAS likely arose from the same signal at the associated locus.

Identification of meQTLs

GWAS was conducted to identify methylation quantitative trait loci (meQTLs) using the significant CpG sites obtained from the epigenome-wide association studies (EWAS) as traits. A total of 2,390 significant associations were identified, corresponding to 2,310 unique CpG sites. To ensure a stringent control for false positive associations, we applied a Bonferroni threshold ($P<1.38\times {10}^{-6}$) for calling significant associations. As a result, we identified 235 associations that surpassed this threshold, involving 52 unique CpGs and 233 unique SNPs (as shown in Table S4).

Among the identified meQTLs, a substantial proportion represented trans-meQTLs, indicating associations between CpG sites and genetic variants located at a distance. Additionally, we observed one cis-meQTL on chromosome SSC12 (Table 5). Specifically, for the cis-meQTL-CpG pair (as illustrated in Figure 4), the SNP (WU_10.2_12_3868856) was found to be approximately 133 Kb away from the CpG site (SSC12:3735523), indicating their close proximity to the FASN gene.

Figure 4. Boxplot of the cis-meQTL (WU_10.2_12_3868856- SSC12:3735523).

The boxes indicate the interquartile range (IQR) of values between the 75% (Q3) and 25% (Q1). The centre lines indicate the median value. The bars below and above each box indicate the data in Q1-1.5 x IQR and Q3 + 1.5 x IQR, respectively. The y-axis shows CpG methylation levels, and x-axis shows SNP genotypes.

Table 5. Summary of meQTL identified by significant CpG sites in this study.

Terms	Value
Number of analysed CpGs	2310
Number of analysed SNPs	36,163
Bonferroni threshold	1.38E-06
Number of meQTL	235
Number of cis-meQTL	1
Number of trans-meQTL	234

Discussion

Based on genome and methylome data, we identified 2390 associations between DNA methylation and 17 fatty acid traits using EWAS. Notably, a substantial proportion of these associations were unique to the EWAS and were not detected by Genome-Wide Association Study (GWAS) methods. Furthermore, by taking methylation of significant CpGs from EWAS as phenotypes, we performed meQTL analysis using GWAS data, resulting in the identification of 235 meQTL-CpG associations. Most of these meQTLs were found to be trans-acting, indicating that genetic variations located at distant genomic regions can influence DNA methylation patterns. Interestingly, we also observed a noteworthy cis-meQTL-CpG pair in close proximity to the FASN gene, suggesting that genetic variations near this gene may impact fatty acid metabolism through the regulation of DNA methylation.

This study provided valuable insights into the associations between DNA methylation and fatty acid traits and elucidated the genetic influence on DNA methylation patterns. Additionally, our findings highlight the importance of considering both EWAS and GWAS approaches to capture the comprehensive landscape of molecular associations.

We identified numerous associations between DNA methylation levels and fatty acid traits (Table 4). It is worth noting that the overlap between EWAS and GWAS associations was relatively low, with only 0.38% of the EWAS associations for fatty acid traits being replicated by GWAS in the same panel of pigs. This limited overlap could be attributed to the challenges in detecting associations with small effects using GWAS, as well as the potential impact of low-density SNP coverage, which may result in the failure to detect certain associations. To further investigate the nature of the overlapping associations, we conducted conditional association studies, revealing that the associations identified by both EWAS and GWAS were likely driven by the same underlying signal at the associated locus.

Previous studies have shown that SNPs near the CpG site may alter methylation levels and result in different DNA methylation patterns (Shoemaker et al. 2010; Zhi et al. 2013). These methylation quantitative trait loci (meQTLs) exert their effects on DNA methylation across the entire genome, thereby potentially impacting phenotypic traits. Uncovering the complex gene network underlying genetic variant-trait associations through meQTL analysis holds significant scientific importance. In our study, we conducted a Genome-Wide Association Study (GWAS) to identify meQTLs by treating the significant CpG sites obtained from the EWAS as phenotypes and the SNPs as predictors. Interestingly, we observed that only 2.25% of the significant CpGs were found to be significantly associated with SNPs. This relatively low overlap could be attributed to our focus on using only the significant CpGs from the EWAS as phenotypes. Notably, previous human meQTL studies have reported associations between genetic variation and 28% of variable CpGs in adipose tissue (Grundberg et al. 2013) and 33% in whole blood (Huan et al. 2019), further supporting the link between DNA methylation levels and genetic variation. Although most (99.57%) of the CpG-SNP associations we identified were trans-acting, we found a cis association on SSC12 near the FASN gene, a gene directly involved in fatty acid metabolism (Zhang et al. 2019). This finding suggests that genetic variation may impact the DNA methylation of the FASN gene, thereby influencing fatty acid traits in pigs.

The growing body of evidence highlights the role of epigenetic modifications, particularly DNA methylation, in contributing to phenotypic variation in livestock species (Triantaphyllopoulos et al. 2016; Wang and Ibeagha-Awemu 2020) and supports the potential application of epigenetic biomarkers, particularly DNA methylation in livestock breeding programs. Moreover, a previous study has reported the association between SNPs and differential DNA methylation (Imgenberg-Kreuz et al. 2018). This provided one possible mechanism that SNPs impact gene expression by altering DNA methylation, thereby suggesting the possible application of epigenetic biomarkers in livestock improvement breeding (Maldonado et al. 2019). Thus, in developing new breeding methods, the relationship between DNA methylation biomarkers identified by EWAS and production traits can be considered to quantify the epigenetic contribution to breeding value prediction. (Wang and Ibeagha-Awemu 2020). However, the application of epigenetic biomarkers in pig breeding still poses challenges. There is a need for comprehensive data on the contribution of epigenetic changes to pig production traits, as well as the development of statistical methods and software to support the quantification of the epigenetic biomarkers’ contribution to phenotypic variation. Therefore, further studies are warranted to better understand the epigenetic mechanisms underlying phenotypic variation in pig production.

In our study, we identified a crucial DNA methylation site and its cis-meQTL on SSC12 near the FASN gene by integrating EWAS and GWAS. This finding holds promise as a potential epigenetic biomarker for future pig breeding programs. However, it is important to acknowledge that our study’s small sample size poses a limitation, which may have influenced our ability to detect significant loci to some extent. Although we have made efforts to increase the sample size to a certain degree, the challenges associated with measuring and recording traits relevant to fatty acids in actual production settings, along with the high cost of Reduced Representation Bisulphite Sequencing (RRBS) for large sample sizes, have constrained the size of our study population.

Authors’ contributions

G.T., Q.S., and K.W. designed the experiments. Q.S., K.W., and S.W. analyzed data. K.W., D.C., X.J., Q.S., and Y.Y. wrote and edited the manuscript. Y.J., X.X., and X.Q. provided the resources. W.X. and A.J. supervised the research. All authors have read and approved the manuscript. Thanks to K.W. for assistance with the experiments and to G.T. for valuable discussion.

Disclosure statement

The authors declare that they have no competing interests.

Data availability statement

The data supporting this study’s findings are available from the corresponding author upon reasonable request.

Additional information

Funding

This study was supported by High-performance Computing Platform of Sichuan Agricultural University. This study was also supported by grants from the Sichuan Science and Technology Program (2020YFN0024), the earmarked fund for the China Agriculture Research System (CARS-35-01A), the National Key R&D Program of China (2018YFD0501204), the National Centre of Technology Innovation for Pigs (NCTIP-XD/B01), the National Natural Science Foundation of China (C170102) and the Sichuan Innovation Team of Pig (sccxtd-2021-08).