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 group, respectively (Table 1). 126,098,724 and 126,782,896 reads were uniquely mapped to the reference chicken genome (Gallus gallus-5.0) in control and inoculated groups (Table 1). The coverage analysis revealed that approximately 81% of the chicken genome was covered by reads at least 1-fold, nearly 77% of genome was covered by more than 5-fold and 55% of genome was covered more than 10-fold (Table 2). In addition, 205,500,619 and 215,922,395 methylated cytosines were detected from control and inoculated group, respectively (Table 3).
The number of methylated cytosines in each type of mCHG, mCHH, and mCpG was counted and the ratio was calculated (Table 3). 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.00% (mCHH/CHH) and 1.00% (mCHG/CHG) identified in the control group. 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. Furthermore, CAG 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 (Fig. 2). The promoter region (the 2 kb bases upstream from the transcription start site (TSS)) had the lower methylation level. TSS had the lowest methylation level. The level of DNA methylation in the first exon was the lowest across all exons, but higher than that in introns. In general, the methylation density in gene body was higher than that in the upstream and downstream of the gene.
Differentially methylated cytosines (DMC) in different genes
DMC was analyzed through MOBAS according to the binomial distribution combined with the 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. The density of miRNAs was higher than other genes. There were 457 miRNAs in the top 1,000 genes with higher 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 the 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%. Also, the ratio on Chr16, Z and W were the lowest ones with 0.07%, 0.01% and 0.002%, respectively. 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. 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 Chr22, 24, 25, 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 a 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. Of 8,946 DMGs, 7,362 genes 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. There were 54.34% (263/484) of methylated genes mapped to immune system process and 50.74% (2241/4417) of methylated genes mapped to metabolic process (Table S3). In term 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 molecular transducer activity, receptor activity, nucleic acid binding transcription factor activity, guanyl-nucleotide exchange factor activity and chemoattractant 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 WNT5B, WNT2 and MYC were involved in Wnt signaling pathway. CCL4, IL-5, IL-6, IL-21 and IL-22 were involved in the pathway of Cytokine-cytokine receptor interaction (Table S4).
Validation of DMGs through Bisulfite Sequencing PCR (BSP)
To validate the sequencing results, six DMGs of HOXA3, HOXD12, DNAH7, NPAT, MTR and ZFHX3 were randomly selected. 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 4). The methylation level detected using BSP method was consistent with that in WGBS results.