DNA methylomic homogeneity and heterogeneity in muscles and testes throughout pig adulthood

The distribution of global CpGs and methylation level

In total, we acquired 2.24Tb WGBS data with an average of 46X genome coverages in 20 samples for BMX and LW, of which 55,685,213 autosomal CpGs (both strands of DNA) were detected. These CpG sites are unevenly distributed across the genome, with the preference of locating at the ends of chromosomes, and abundantly clustered in local locations (Figure 1A).

Integrating genome annotation information, we detected the proportion of these CpG sites distributed in introns, exons, and intergenic regions are 51.9%, 9.11%, and 38.98%, respectively (Figure 1B). Further analysis uncovered that 9.62% of 55,685,213 CpG sites are enriched in CpG islands (CGIs), of which 11.64% CpGs are located in CGI shores (Figure 1C).

We retained about 35M CpGs whose coverage ≥10 and ≤300 in half of samples to investigate the global methylation level. The global methylation levels formed a skew bimodal distribution with mean 70.5% (Figure 1D). In LD, the global methylation level was higher in LW than in BMX. Comparing results between tissues, the methylation level in testis was higher than in LD for BMX (Figure 1E). Based on the annotated genes and CGIs in pig genome, we also showed the distribution of methylation levels over ten regions (Figure 1F). Lower mean methylation levels (~50%) were found in 5’ UTRs and CGIs, and methylation levels in the two regions had a larger inter-quartile (~0.8). In addition, overlap of 5’ UTRs and CGIs (5’ UTR-CGI) had an extremely low methylation level (mean 7.65%) (Figure 1F and Supplementary Table 1). Focusing on hypomethylation and hypermethylation CpGs, we explored the distribution of 1,147,532 hypomethylated CpGs (methylation level was lower than 10%) and 5,743,199 hypermethylated CpGs (higher than 90%). More than half of hypomethylated CpGs are localized in CGIs or CGI shores, and the hypomethylated sites are also enriched in 5’ UTRs (Fisher’ exact test odds=2.80, P=2.2e-16). Our results revealed hypermethylated global genome and local hypomethylation in genome and further discovered highly variable methylation level and hypomethylation in 5’ UTRs and CGIs.

To further investigate relation between gene expression and methylation, we focused on the correlation between methylation level in the promoter regions and gene expression as most of CGIs are located in gene promoters in mammalian. We herein calculated the mean methylation level in 1000 bp promoter regions (800 bp upstream of TSS to 200 bp downstream of TSS). We retained 7558 genes’ promoters which have at least ten CpGs in our filtered data. About 6.8% (655 out of 7558) mean methylation level of promoters were strongly correlated with genes’ expression (Spearman correlation, P<0.01, mean |r|=0.66). The results suggested that gene promoter methylation plays an importantly regulatory role in gene expression in the adult tissues.

Differentially methylated CpGs between tissues

To reveal tissues or organs specific DMCpGs relating to structural and functional elements in adulthood pigs, we detected the methylomic differences between LDs and testes from young to elder stages. The joint results of three differential analyses under BGLMM framework identified ~5.3 million DMCpGs (9.52% out of captured autosomal CpGs, P<0.01, Figure 1A) which belonged to 1,466,027 discrete DMRs following Matthew D.S et al. [15]. Focusing on the DMRs in promotors, we uncovered 110582 DMRs positing within 1^st exon or 1Kb upstream of genes and rigorously screened out 4978 DMRs (53627 DMCpGs) which included at least 5 synclastic DMCpGs. Correlation analysis showed 7361 DMCpGs and 5197 DMCpGs were negatively and positively correlated with genes’ expression, respectively (mean correlation coefficient was -0.59, P<0.05 and 0.57, P<0.05). The methylation status of DMRs in adult tissues may have positive and negative effects on the expression levels of 8% and 11% of genes (Supplementary Figure 3). DMRs of 1362 genes were detected to be significantly differentially expressed between two tissues (Wilcoxon signed rank test, P<0.01, Supplementary Table 2), among which DMRs locating in genes’ promoters displayed higher differential expression tendency (Fisher’s exact test, odds ratio=1.27, P=2.82E-12). There were no significant differences in expression of the remaining 3616 genes. These results implied that a part of DMRs continuously regulated the tissue-related function through all adulthood, and others might be of spatiotemporal characteristics or other unknown regulatory mechanisms.

For the known functional and significantly differentially expressed 873 genes, we performed Gene ontology (GO) analysis and found that 549 genes were enriched into 101 GO terms (P<0.05) in three categories including biological process (65 terms), cellular component (19 terms) and molecular function (17 terms). Muscle tissue-related function term includes muscle contraction, muscle filament sliding, myofibril assembly and skeletal muscle cell differentiation. Testis specific functional ontology includes meiotic DNA double-strand break formation, sister chromatid cohesion, spermatid development and acrosomal membrane (Supplementary Table 3). In addition, 70 genes were enriched into a biological process of positive regulation of transcription from RNA polymerase II promoter (GO:0006936). The results showed that DMCpGs were associated with tissue-specific cell differentiation and gene expression.

The mechanisms of DNA methylation in regulating gene expression are very complex and diverse. Herein we report three different DMRs in promotor whose methylation levels are negatively correlated with gene expression (Figure 2A–2C), including DMR co-localized with CGI (SPACA1), differential DNA methylation valley nearby the transcriptional start site (MEI1), and DMR co-localized with conservative sequence (TNNT3). Besides, we also detected a certain number of positive correlations between methylation in upstream proximal DMRs (TCEA3) and its expression (Figure 2D). It supports the current research results that methylated regulatory elements activate gene expression [10, 11].

Normally, gene expression is regulated by gene promotor methylation. Interestingly, DMR on chromosome 2 (chr2:1301111-1301408) correlated with expression was located in TNNT3 rather than in its promotor region (Supplementary Figure 2). This region is probably not the transcription initiation site in muscles based on RNA-seq data (Figure 2C). Conservative analysis for 21 eutherian mammals revealed an 83bp constrained EPO-low-coverage element (chr2:1301181-1301263) was co-localized in the DMR, and the methylation levels of CpGs between constrained region and the ~100bp flanking region were significantly different in two tissues. Thus, the 83nt core sequence has been browsed to AnimalTFDB3.0 (http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/), and 64 pieces of information was retrieved (P<10^-4). However only five expressing predicted transcription factor genes (including NFATC2, NFATC4, RBPJ, RB1 and TEAD4) were identified in the two tissues (Figure 2E). The results showed the constrained element is an alternative promotor or tissue-specific enhancer in TNNT3 of pig.

Dynamic CpGs with age

In order to clarify the physiological tendency in LD and testis under three age status, the sections of muscles and testes were stained with ATPase and HE, respectively (Figure 3A, 3B). Histological sections of LDs showed reducing of slow muscle fiber ratio and hypertrophic multinucleated muscle cells with age (Figure 3B). Histological sections of testes showed a decrease in the number of germ cells and Sertoli cells, and a shrinkage in seminiferous tubules with aging (Figure 3A). the 9Y testes showed typical histological aging characteristics [28].

We herein compared the methylation level among ages to see that if the distinct phenotype was regulated by methylation. Our methylomic data showed no significant difference in methylomic entropy and global methylation level across age groups (ANOVA, P>0.05), but weak decreasing tendency of methylation level was observed (methylome level in 1Y, 4Y, and 9Y age groups were 72.94%, 70.95%, 69.72%, respectively).

Focusing on some specific CpGs, we detected 9821 and 9759 dynamic CpGs with age using BGLMM in LDs and testes respectively. The dynamic CpGs are unevenly distributed on the genome, of which 766 (~7.8%) have been found in both tissues. The CpGs varied consistently in both tissues, which suggested methylation hotspot with ageing (Figure 1A). The dynamic CpGs showed more hypomethylation with age (57%) than hypermethylation (43%), and 69% of shared dynamic CpGs were hypomethylated with age. The Short Time-series Expression Miner (STEM) analysis revealed several differentially dynamic patterns in CpG methylation with age (Supplementary Figure 4).

Using methylation difference threshold algorithm (minimum methylation difference ≥0.3), we identified 147159 and 294439 dynamic CpGs varying with age in LD and testis, respectively (Figure 3C, 3D). We validated 56% dynamic CpGs identified by BGLMM algorithm. A total of 43,655 dynamic CpGs were shared in the two tissues and more hypermethylation with age than hypomethylation was showed. The lack of regular MDCpGs supports the previous hypothesis of stochastic epigenetic drift with aging. However, we captured a gene set including 93 genes which harbored at least 50 intragenic DMCpGs in both tissues. The shared genes enriched into 14 biology process terms which refer to neuron migration and development, signal transduction, and cell adhesion (Figure 3E). Among them, at least 16 homologous genes (including AKT3, ARNT, CDC42BPA, CDH4, COL11A1, CTNNBL1, MYO1D, PRKCE, PSD3, PTPN4, RAD51B, RGS7, SLC24A3, TGFB2, TMEM117, UNC5D) are listed in the aging gene database (http://genomics.senescence.info/genes/).

Differentially methylated CpGs in LD between breeds

Focusing initially on our muscle samples of two pig breeds, we captured 1147591 dynamic CpGs, of which 63.31% are distributed within 25959 genes. We detected 7876 co-directional methylated CpGs in 605 genes which were significantly differential expressed (Wilcoxon signed rank test, P<0.01) in LDs between two breeds. The HOXC gene cluster (ssc12:24764000-2489000) exhibited rich differential methylation sites and 38.77% of which were located in CGIs. Further analysis of RNA-seq data indicated 6 out of 9 HOXC family genes were differentially expressed (Figure 4), the expression levels of HOXC4-6, HOXC8, HOXC9 in LD of LW were higher than that of BMX, and the HOXC11 was in opposite in LD of two breeds (Figure 4). This results suggested that intragenic and intergenic DMRs may affect expression of the host gene. We focused on HOXC8 which was highly expressed in LD of LW. Muscles may require HOXC8 protein for full activation of muscle-specific gene expression [29]. Up regulation of HOXC8 in LD suggests its role in lineage development of muscle satellite cell (MSC) into trunk muscles [30].

Other two genes, ZIC1 and BMP5, were highly expressed in LD of LW. BMP signaling regulates MSC-dependent postnatal muscle growth [31], and ZIC1 is a marker of adipose tissue browning in humans, which is defined as the conversion of white fat into brown fat [32]. Compared with LW pigs, BMX pigs have lower eye muscle area (13.63 cm² vs 51.47 cm²) and higher intramuscular fat content (IMF, 2.93% vs 1.89%) [33, 34]. The results suggested a set of highly specific-expressed genes with hypomethylation in LD of LW are beneficial for LD hypertrophy and maintenance and are not conducive to the deposition of IMF.

DMRs between two pedigrees

The breed-specific methylomic pattern suggested that partial DNA methylation characteristics were stably transmitted within single breed. To verify intergenerational inheritance, we further compared the methylation levels between two pedigrees, and a total of 64 and 17 inheritable DMRs were identified in testes and muscles, respectively (Adjusted false positive <0.01, Methods, Supplementary Table 4). The parental methylation patterns in genome regions are transmitted to the offspring, strongly supporting the heritability of partial methylated or demethylated sites in adult tissues (Figure 5). We carefully analyzed an about 130bp DMR on SSC9 in testes which harbored 14 DMCpGs. After correcting the effect of sequencing depth on methylation level by the posterior probability of beta distribution [14], the methylation pattern differences between pedigrees are more than 0.3, suggesting intergenerational or transgenerational epigenetic inheritance.

Dynamic CpGs between genders

We found that 155211 of CpG sites (567289 DMRs) in these experiments are strongly differentially methylated between genders (Minimum methylation difference ≥0.3), of which 69 are located in gene promotor DMR and 3419 are intragenic. Differential expression tests for DMR genes suggest no significant or weak differences (Wilcoxon signed rank test, P>0.05) in expression level whether they are in the promoter or intragenic regions.

Animals and sampling

All procedures involving animals are in compliance with guidelines for the care and use of experimental animals established by the Ministry of Agriculture of China, and the trial was approved by the Animal ethics committee at Jiangxi Agricultural University (No. JXAULL-2016001).

This study involved two pig breeds with dramatic genetic differences, one is western commercial purebred of Large White (LW) and the other is Chinese indigenous purebred of Bamaxing (BMX). Four LW pigs (about seven months, two virgin female and two male pigs) and eight male BMX pigs (including three nine-year-old (9Y), three four-year-old (4Y), and two one-year-old (1Y) pigs). They were raised at the farm and were fed on a similar diet, and they could access water and food ad libitum under a standardized feeding and management regimen. The pedigree information of BMX is showed in Supplementary Figure 1. After slaughting within postmortem 45 min, about 1 gram (g) portion of Longissimus dorsi (LD) muscle at the 1^st-2^nd lumbar vertebra on the left side of the carcass and about 1g testis at central position of left testis were sampled and rapidly stored in −80° C refrigerator for DNA extraction.

Whole genome bisulfite sequencing

A total of 2.26 Tb (average 113.20 Gb) WGBS data were generated using Hiseq X10 platform for twelve skeletal muscles and eight testes. The average mapping coverage and depth rate were 2.43 Gb and 46.49X per sample, respectively. The clean reads were aligned to the Sus scrofa 11.1 using bwa-meth-master, and an average of 96.92% reads were uniquely aligned to the reference genome. MethylDackel0.3.0 was used to summary methylation state and a total of 55,685,213 autosomal CpGs have been identified in at least one of twenty samples. These CpGs were annotated using Sus_scrofa.Sscrofa11.1.98.gtf which was downloaded from Ensembl (http://asia.ensembl.org/index.html).

RNAseq

An average of 10G RNA-seq clean reads using NovaSeq5000 platform from 20 samples were mapped to Sus scrofa 11.1 with hisat2.1.0. The mapped reads were further analyzed by StringTie1.3.6 and subread1.6.5, and the expression levels for each transcript were quantified as fragments per kilobase of transcript per million mapped reads (FPKM).

Differential analysis

CpGs with the coverage of less than 10 and greater than 300, and CpGs with individual missing rate less than 50% were filtered out. Three methods including MACAU [44], MALAX [45], and PQLseq [46] have been introduced for differential analysis. These methods adopted a binomial generalized linear mixed model (BGLMM) for count-based bisulfite sequencing data, but they used different algorithms including markov chain mote carlo (MCMC) sampling-based strategy [44], Laplace approximation [45], and penalized quasi-likelihood approach [46], respectively. The model treats age as a covariate and considers the covariance among individuals due to individual relatedness (Supplementary Table 5). The significant threshold was set as 0.01, the intersection result of three methods was used as DMCpGs. DMR was defined as combined sites within 500 bp of one another [15].

We also adopted a previously published method [14] which is based on the beta distribution to model single CpG methylation levels and exhibited a minimum methylation level difference of ≥0.3 at a significant level of P≤0.01. We estimated the posterior distribution of methylation levels in both genders from LW and three age groups from BMX [14], and we tested the minimum methylation level differences.

Functional annotation

Pig methylation sites were annotated using Sus_scrofa.Sscrofa11.1.98.gtf, and the annotation of CGIs of Sus scrofa 11.1 reference sequence was downloaded from the University of California Santa Cruz (UCSC) Genome Browser. It was annotated following a criterion of segment length >200 bp, CG content >50%, and observed/expected ratio of CpG sites >0.6 [47]. CGI-shores are defined as the ~2 kb regions near CGI, and CpG shelves are further 2kb extension regions of CGI-shore. Differentially methylated gene lists were submitted to online tool DAVID (https://david.ncifcrf.gov/). We utilized homo sapiens set as “background” and selected GO terms and KEGG terms based on the statistical significant level (P<0.05).

Inheritable methylation regions

The heritable region was defined as a 500 bp DMR in which at least 5 synclastic DMCpGs were detected for all individuals between pedigrees. In statistics, it is a very strict threshold (P_fp value after Bonferroni correction is 0.01) to exclude random false positive. P_fp value was calculated as follows:

The number of DMCpGs (N) within segments (L bp) follows a binomial distribution Bin(N, ECG, p/2).

$$P(X=N)=C_{ECG}^N{\left(\frac{p}{2}\right)}^N{\left(1-\frac{p}{2}\right)}^{ECG-N}$$

The expected number of CpGs (ECG):

$$ECG=\frac{NumberofC\times NumberofG}{L}$$

The false positive probability by Bonferroni correction:

$$P_{fp}(X\ge N)=2(1-P(X<N))\frac{GS}{L}$$

The ratio of DMCpG (p) is less than 0.003 between pedigrees in our data; The expected number of CpGs (ECG) is about 31 in a 500bp (L) random sequence; Pig genome size (GS) is about 2.5G. Under the ECG, the Bonferroni correction P value is:

$$P_{fp}(X\ge 5)=0.012$$