Analysis of genome-wide DNA methylation data
One genomic DNA library was constructed in each of control and inoculated groups. There were 185,362,463 and 180,530,750 clean reads obtained from control and inoculated library, respectively (Table 1). Of those clean reads, 126,098,724 and 126,782,896 were uniquely mapped to the reference chicken genome (Gallus gallus-5.0) (Table 1). In addition, 205,500,619 and 215,922,395 methylated cytosines were detected from control and inoculated group, respectively (Table 2).
Table 1
Data generated by whole genome bisulfite sequencing
Sample | Clean reads | Clean Base | Unique mapped reads | Mapped (%) | Conversion rate (%) | GC(%) |
Control | 185,362,463 | 55,494,074,424 | 126,098,724 | 68.03 | 99.04 | 23.59 |
Inoculated | 180,530,750 | 54,074,103,970 | 126,782,896 | 70.23 | 99.25 | 23.12 |
Table 2
Number and ratio of different types of methylated sites
Group | mCHG | mCHH | mCpG | Total mC |
Control | 14,037,071(1.00%) | 37,592,584(1.00%) | 153,870,964(55.20%) | 205,500,619 |
Inoculated | 13,045,463 (0.90%) | 35,035,989(0.82%) | 167,840,943(55.20%) | 215,922,395 |
*H = A/T/G |
The number of methylated cytosines in each type of mCHG, mCHH, and mCpG was counted and the ratio was calculated. There were 153,870,946 mCpG sites, 37,592,584 mCHH sites and 14,037,071 mCHG sites accounting for 55.20% (mCpG/CpG), 1% (mCHH/CHH) and 1% (mCHG/CHG) identified in the control group (Table 2). 167,840,943 mCpG sites, 35,035,989 mCHH sites and 13,045,463 mCHG sites accounting for 55.20% (mCpG/CpG), 0.82% (mCHH/CHH) and 0.90% (mCHG/CHG) were identified in the inoculated group (Table 2). Furthermore, CAG context was dominant in mCHG type. CAH and CHT were preferred in the mCHH type (Fig. 1).
DNA methylation in different gene regions
To better understand methylation pattern in the genome, the methylation in different gene regions was analyzed. Promoter region (the 2 kb bases upstream from the transcription start site (TSS)) had the lower methylation level (Fig. 2). TSS had the lowest methylation level (Fig. 2). The level of DNA methylation in the first exon was the lowest across all exons, but higher than that in introns of gene body (Fig. 2). In general, the methylation density in gene body was higher than that in upstream and downstream of gene.
Characterization of differentially methylated cytosines (DMC) in different genes
DMC was analyzed through MOBAS according to the binomial distribution combined with bayesian algorithm. The distribution of DMCs in different chromosomes was shown in Fig. 3. The density of DMC located in Chr1-3 was lower than that in Chr5-30. And the density of DMC in Chr W and Z was the lowest (Fig. 3). The density of miRNAs was higher than other genes. There were 457 miRNAs in the top 1,000 genes with high DMC density, 324 miRNAs in the top 500 genes. Gga-miR-7466, gga-mir-1713, gga-mir-1699, gga-mir-7467, gga-mir-6616 had the highest methylation density (Table S1). The HOX gene family was widely methylated and mainly distributed in Chr2 and 7 (Table S2).
Identification of differentially methylated region (DMR) and differentially methylated genes (DMG)
There were 82.5% DMR located in distal intergenic region, only 0.02% in the 1st intron (Fig. 4). The DMR coverage ratio on each chromosome was calculated. The coverage ratio on Chr1, 2, 3, 16, 25, 31, 32, 33, Z and W was less than 0.5%, ratio on Chr4, 5, 6, 7, 8, 11, 22, 30, 31 was between 0.5% and 0.8% and the ratio on Chr9-28 excluding for Chr11, 16, 22 and 25 was more than 0.8% (Fig. 5). In addition, the ratio on Chr16, Z and W were the lowest ones with 0.07%, 0.01% and 0.002%, respectively (Fig. 5). The ratio on Chr23, 24 and 26 were the highest ones with 1.10%, 1.18% and 1.22%, respectively (Fig. 5).
There were 8,946 differentially methylated genes identified including 3,639 hypo-methylated genes and 5,307 hyper-methylated genes in the inoculated group compared with the control group. Differentially methylated genes distributed variously across all chromosomes (Fig. 6). There were more than 1,000 differentially methylated genes in Chr1, 500-1,000 differentially methylated genes in Chr2, 3, 4 and 5, 100–500 differentially methylated genes in the Chr6-28 excluding Chr16, 22,24 and 25, 10–100 differentially methylated genes in Chr16, 22, 24, 25, 32, 33 and Z. Number of differentially methylated genes in Chr16, 32 and W was less than 10. More hyper-methylated genes were identified than hypo-methylated genes on all chromosomes except for Chr16 and W.
COG function classification of differentially methylated genes
The COG (Clusters of Orthologous Groups) function classification results showed that the DMGs were mainly associated with seven categories: general function prediction only, signal transduction mechanisms, transcription, replication, recombination and repair, posttranslational modification, protein turnover, chaperones, amino acid transport and metabolism and inorganic ion transport and metabolism with percentage of 39.56, 16.12, 15.48, 14.51, 7.65, 6.18 and 5.79%, respectively (Fig. 7).
Functional annotation of differentially methylated genes
To understand the function of those differentially methylated genes, Gene Ontology (GO) and KEGG pathway enrichment were analyzed. Totally, 7,362 of 8,946 DMGs were annotated. The results of BP (biological processes), MF (molecular functions), and CC (cellular components) were shown in Fig. 8. For the BP, the DMGs were mainly associated with immune system process, metabolic process, reproductive process, signaling, multicellular organismal process, developmental process, hormone secretion, rhythmic process, response to stimulus, biological regulation and cell aggregation (Fig. 8). There were 54% of genes mapped to immune function term methylated, and 50.74% of genes mapped to metabolic process term methylated (Table S3). In terms of the CC, the DMGs were mainly located in the extracellular region, cell, nucleoid, organelle part, virion part and membrane part. For the MF, the DMGs were associated with translation regulator activity, nutrient reservoir activity and protein-binding transcription factor activity.
There were 16 KEGG pathways associated with DMGs significantly enriched (P < 0.05). The enriched pathways were roughly grouped into three groups: 1) immune-related pathways including Cytokine-cytokine receptor interaction, TGF-beta signaling pathway, FoxO signaling pathway, MAPK signaling pathway and Wnt signaling pathway; 2) metabolism-related pathways including mTOR signaling pathway, other types of O-glycan biosynthesis, inositol phosphate metabolism, Glycosphigolipid biosynthesis-lacto and neolacto series, Alanine, aspartate and glutamate metabolism, Mucin type O-Glycan biosynthesis, and Glycosaminoglycan biosynthesis-heparan sulfate/heparin, Melanogenesis; 3) others including Progesterone-mediated oocyte maturation, adherents junction, vascular smooth muscle contraction, neuroactive ligand-receptor interaction (Fig. 9). There were 92, 80 and 98 DMGs associated with Wnt signaling pathway, FoxO signaling pathway, and Cytokine-cytokine receptor interaction, respectively. Genes in Wnt family like Wnt8A, Wnt8B were involved in Wnt signaling pathway. CCL4, IL-5, IL-6, IL-21, IL-22 were involved in Cytokine-cytokine receptor interaction (Table S4).
Validation of DMGs through Bisulfite Sequencing PCR (BSP)
To validate the sequencing results, six DMGs were randomly selected. The four loci in HOXD12 gene (Fig. 10) and two loci in MTR gene were hyper-methylated, two loci in ZFHX3 were hypo-methylated in the inoculated compared with the control group (Table 3). The loci in Chr2:32659632, Chr7:9749766, Chr1:180358709 and Chr1:180358843 were hypo-methylated, while loci in Chr2:32659658, Chr7:9750310 and Chr1:180358844 were hyper-methylated (Table 3). These validations were basically consistent with the WGBS data.
Table 3
The methylation level of validated genes between GWBS and BSP method
Gene name | Location | Meth_direction in GWBS | BSP |
Inoculated group | Control group | Difference/ meth_direction |
HOXA3 | Chr2:32659632 | Hypo- | 0.4 | 0.6 | -0.2 |
Chr2:32659658 | Hyper- | 0.7 | 0.2 | 0.5 |
HOXD12 | Chr7:16384244 | Hyper- | 1 | 0.2 | 0.7 |
Chr7:16384255 | Hyper- | 0.7 | 0.4 | 0.3 |
Chr7:16384275 | Hyper- | 0.5 | 0 | 0.5 |
Chr7:16384301 | Hyper- | 0.8 | 0 | 0.8 |
DNAH7 | Chr7:9749766 | Hypo- | 0.4 | 0.8 | -0.4 |
Chr7:9750310 | Hyper- | 0.8 | 0.2 | 0.6 |
NPAT | Chr1:180358709 | Hypo- | 1 | 0.7 | 0.3 |
Chr1:180358843 | Hypo- | 0.6 | 0.9 | -0.3 |
Chr1:180358844 | Hyper- | 0.9 | 0.5 | 0.4 |
MTR | Chr3:37701958 | Hyper- | 0.7 | 0.2 | 0.5 |
Chr3:37702200 | Hyper- | 0.6 | 0.1 | 0.5 |
Chr3:37702201 | Hyper- | 0.7 | 0 | 0.7 |
ZFHX3 | Chr11:19406308 | Hypo- | 0.2 | 0.5 | -0.3 |
Chr11:19406309 | Hypo- | 0.3 | 0.7 | -0.4 |