Cecum methylome response to Salmonella enterica serovar Enteritidis infection in chicken through whole genome bisulte sequencing

Background: Salmonella enterica serovar Enteritidis (SE) is one of the pathogenic bacteria, which affects poultry production and poses severe threat to public health. Chicken meat and egg are the main source of SE. DNA methylation, an important epigenetic modication, involves in regulatory processes including gene expression, chromatin structure and genomic imprinting. To understand the methylation regulation in response to SE inoculation in chicken, the genome-wide DNA methylation prole following SE inoculation was analyzed through whole genome bisulte sequencing in the current study. Results: There were 185,362,463 clean reads and 126,098,724 unique reads in the control group, and 180,530,750 clean Reads, 126,782,896 unique reads in the inoculated group. We found that the methylation density in gene body was higher than that in the upstream and downstream regions of gene. There were 8,946 differentially methylated genes (3,639 hypo-methylated genes, 5,307 hyper-methylated genes) obtained between inoculated and control groups. Methylated genes were mainly enriched in immune-related Gene Ontology (GO) terms and metabolic process terms. Cytokine-cytokine receptor interaction, TGF-beta signaling pathway, FoxO signaling pathway, Wnt signaling pathway and several metabolism-related pathways were signicantly enriched. The density of differentially methylated cytosines in miRNAs was the highest. HOX genes were widely methylated and mainly distributed in Chr2 and 7. Conclusions: We rstly analyzed the genome-wide DNA methylation in the response to SE inoculation in chicken. SE inoculation promoted the DNA methylation in chicken cecum and caused methylation alteration in immune- and metabolic- related genes. Wnt signal pathway, miRNAs and HOX gene family may play a crucial role in the methylation regulation of SE infection in chicken. The ndings herein will deepen the understanding of epigenetic regulation in the response to SE inoculation in chicken.

Background Salmonellosis, mainly caused by Salmonella, is one of the most frequently infectious foodborne diseases around the world. Salmonella enterica serovar Enteritidis (SE) is one of pathogenic bacteria, which causes signi cant economic losses on poultry production and put severe threat on human health [1][2][3] through poultry and egg [4][5][6]. It is estimated that Salmonella enteritidis causes 1.3 million cases of gastroenteritis and more than 350 died each year in the United States [7]. In Europe, 96,039 salmonellosis cases were reported in 2016 [8]. Effective prevention of SE infection in poultry production has aroused researchers' attention.
It has been reported that SE infection affects gene expression in chicken [9][10][11]. The expression level of Toll-like receptor 2 (TLR2) and TLR4 increased signi cantly in chicken spleen infected with Salmonella enterica serovar Enteritidis [12]. The immune response was suppressed when inoculated with SE at the onset of laying in White Leghorn layer [13]. Moreover, miRNAs contribute to the response to SE infection at the onset of egg laying through regulating the homeostasis between metabolism and immunity [14].
Epigenetic modi cation plays a critical role in genomic stability and gene expression. DNA methylation, one of major epigenetic modi cations found in most of the eukaryotes, could participate in many biological processes, like gene expression regulation, development and disease-resistance [15][16][17].
DNA methylation in gene body interferes with transcript elongation [18][19][20], and DNA methylation of the rst exon is tightly linked to transcriptional silencing [21,22]. Genome-wide DNA methylation mapping uncovers epigenetic modi cation changed with animal development, evolution and environmental adaptation [15,[23][24][25]. Genome-wide DNA methylation pro les have been reported in many species such as human [26], bovine [15], soybean [27], rat [28], rice [29] and chicken [30]. DNA methylation is implicated in human disease and immune system [31] through regulating the gene expression [32,33]. Interleukin 6 (IL-6) promoter in peripheral blood mononuclear cells is hypomethylated in Rheumatoid Arthritis patients [34]. IL-10 promoter was differentially methylated in arthritic cells [35]. The promoter of IL-1R2 was hypomethylated in the rheumatoid arthritis patients [36]. It has been reported that an arthritis-speci c DNA hypomethylation event could regulate autoimmune arthritis by interfering with an anti-in ammatory pathway [37]. The hypomethylation of ZBT38 in B cell represses expression of IL-1R2 gene and lead to compromised anti-in ammatory [37]. Furthermore, aberrant DNA methylation has been associated with several immune de ciencies and autoimmune disorders [38]. Multiple functional DNA methylated loci have been identi ed to play important roles in the regulation of gene expression involved in the in ammatory response and tissue injury after APEC infection in chicken [39]. Our previous results showed that SE infection repressed overall genomic DNA methylation level in chicken through Methylated DNA quanti cation kit [40]. However, the genome methylation pro le changes following SE inoculation in chicken is still unclear. Whole genome bisul te sequencing (WGBS) [41] has been widely used in studies of DNA methylation associated with growth [42], development [43] and disease [44], allowing for unbiased genome wide DNA methylation pro ling. In the current study, the whole genome-wide DNA methylation pro le following SE inoculation in chicken cecum was investigated through WGBS.

Results
Analysis of genome-wide DNA methylation data   (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 rst 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).

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.

Validation of DMGs through Bisul te 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). These validations were basically consistent with the WGBS data.

Discussion
Salmonella enterica serovar Enteritidis poses sever threat on public health and has been widely studied in different levels. DNA methylation plays a crucial part in the regulatory processes including gene expression, chromatin structure, genomic imprinting, transposon silencing, X-chromosome inactivation, disease response and individual development [45][46][47]. The DNA methylation was identi ed to be associated with SE infection in chicken [48]. The cecum is the primary fermentation site of Salmonella.
We focused on the genome methylation variation of chicken cecum following SE inoculation three days and aimed to reveal the relationship between DNA methylation and SE inoculation.
The DNA methylation mainly occurred in CG context but rarely in non-CG (CHG and CHH) context on all chromosomes or genome functional regions of chicken, which was consistent with the DNA methylation pro le in mammals [49]. It was reported that hypo-methylation in promoter and hyper-methylation in gene body could affect transcription positively, whereas methylation in the promoter generally induces repression of transcription [50]. In the current study, DNA methylation level in gene body was higher than that in gene initiation and gene termination, which was consistent with the previous results [51,52]. The collectively results suggests the methylation of the gene body regions may play an important role in regulating gene expression.
It has been reported that hyperhomocysteinaemia, dyslipidemia, in ammation and oxidative stress could promote aberrant DNA methylation [53][54][55]. We identi ed the number of hyper-methylated genes in inoculated group was higher than that of hypo-methylated genes across chromosomes rather than Chr16 and W. It suggested that SE incubation promoted DNA methylation in chicken and distributed unevenly across Chromosomes, which was inconsistent with our previous study [40,56] probably due to different method or different host genetic background. Previous studies in human [17] and plant [57] showed that individuals with a higher DNA methylation level in some special genes are susceptible to diseases or bacterial infection. It has been reported that the methylated genes following Marek's disease infection in chicken are associated with response to stimulus, cell adhesion, and immune system process function [58]. Previous studies also veri ed the immune and metabolism-related process were related with abiotic stress [59,60] and various disease [14,39,61,62]. In our study, DMGs were mainly involved in metabolism process and immune system process, which was consistent with previous studies. It has been suggested that mechanism of disease resistance determined by DNA methylation and methylation variation could be both a cause and consequence of viral infection [63]. Variation in the methylation state of in ammatory bowel disease (IBD)-associated genes could alter gene expression, potentially contributing to disease onset and progression [64]. The results collectively suggest that DNA methylation may regulate host immune response via modulating expression of certain immune-related genes in response to SE infection.
Rhythmic process plays a dominant role in determining overall health and physiological homeostasis and could facilitate the organism survive well in different circumstances. There was a direct molecular link existed between circadian dysregulation and diseases [65]. It has been revealed that circadian rhythm related genes play a critical role in host response to C. jejuni colonization in chicken [62]. The circadian rhythm associated genes were not enriched following SE infection [13]. While, the circadian rhythmassociated circRNAs and miRNAs were signi cantly triggered by SE inoculation [14,66]. The methylation variation of rhythmic process related genes was triggered in chicken following SE inoculation. The rhythm related process would regulate SE inoculation in different levels and the different response of circadian rhythm between C. jejuni inoculation and SE inoculation need to be further studied. And immune system may correlate with metabolic system involved in regulating SE inoculation. Interestingly, the differentially methylated genes were found to be correlated with immune process, metabolic process and rhythmic process. Several studies also found circadian clocks, metabolism and immune process are inextricably intertwined [67,68]. There may exist an interaction among the three processes in chicken inoculated with SE, however, the detailed mechanisms need to be further studied.
Wnt signaling pathway is involved in the development, cell differentiation and disease pathophysiology [69]. Wnt signaling pathway was signi cantly phosphorylated in SE infected cecum [70]. He et al found that the porcine WNT10B gene may be a key candidate gene for porcine backfat thickness [71]. Liu reported Salmonella activated the Wnt/b-catenin signaling pathway to regulate stem cells [53]. Wnt ligands have been reported to regulate multiple aspects of intestinal pathophysiology. Wnt signaling and mTOR pathways mediate the impaired epithelial autophagy [72]. The enriched Wnt signaling pathway and mTOR signaling pathway may indicate SE inoculation would affect the chicken development through DNA methylation of genes in those signaling pathways.
HOX genes, a conserved gene family, plays a crucial role in embryonic development and are involved in the reproduction and development of cells [73]. The immune function of HOX genes has been reported in tumor, cancer [74]. HOX genes methylation was identi ed to be associated with hereditary breast cancer [75]. A subset of progressively more posterior HOX genes, which are collinearly activated in vertebral precursors, repress Wnt activity with increasing strength [76]. For example, the distinguishing feature of HOXA5 in acquiring and establishing the epigenetic system was mediated through DNA methylation induced gene regulation and also provided a cue for mouse evolutional initiation [73]. Aubin et al found the loss of HOXA5 gene function in mice could perturb intestinal maturation [77]. HOXC10 was identi ed to play an important role in regulation network of growth and reproduction in Jinghai Yellow chicken [78]. Raval et al reported on the methylation of DAPK1 was repressed and interacted with HOXB7 in both sporadic and familial chronic lymphocytic leukemia (CCL) [79]. Tong et al detected HOXB7 gene in CLL is hypermethylated [80]. HOX gene family is a big family, the DMC was detected in 37 genes in HOX gene family, 8 of those genes had higher density than 0.005 and distributed on both Chr2 and Chr7. We speculated SE inoculation could affect chicken growth and reproduction through regulating the methylation of HOX genes.
MiRNA and DNA methylation are two categories of epigenetic modi cations which regulate gene expression. DNA methylation and histone modi cations could regulate miRNA expression during tumorigenesis [81]. In breast cancer cell lines, almost one-third of miRNA promoters are hypermethylated relative to normal mammary cells [82]. Mir-129-2 and mir-663a were highly methylated in human urothelial carcinoma [83]. Methylation of miR-9 is a biomarker for oral and oropharyngeal squamous cell carcinomas [84]. Methylated miR-17-5p was detected in all pancreatic cancer patient samples and is suitable to be a biomarker [85]. Tang found that methylation-sensitive mir-345 play an important role of antineoplastic as a growth inhibitor in the development of colorectal cancer [86]. Li et al found that the methylation level of miRNA promoters was high in chicken infected with Marek's disease [61], which was consistent with our study. The hypermethylation in upstream region of gga-miR-130b-3p gene contributed to its repressed expression in tumorous tissues [87]. In the current study, miRNA has the most density DMCs among all genes. MiR-138 is down-regulated in human colorectal cancer tissues and cell lines [88].
The ndings indicated that miRNAs may be sensitive to DNA methylation responding to SE inoculation in chicken. The tRNA modi cations affect all aspects of tRNA biology including decoding and charging e ciency and delity, in vivo stability, and intracellular localization [89]. Methylation level of tRNA in Escherichia coli and Salmonella sensitize these bacteria to antibiotics [90]. The DMC of 6 tRNAs was higher than 0.10 in the current study. We respected those methylated miRNAs and tRNAs may be referred as a maker in detecting the SE infection in chicken.

Conclusions
In conclusion, the cecum genome-wide methylation pro les of Jining Bairi chicken following SE inoculation were studied through the whole genome bisul te sequencing. SE inoculation promoted the genome-wide methylation level in Jining Bairi chicken. SE inoculation would trigger the aberrant methylation of genes in Wnt signaling pathway, mTOR signaling pathway, immune and metabolism related functional terms. Methylation of miRNA and HOX gene family may play roles in the epigenetic regulation responding to SE inoculation. Our results would pave the foundation for understanding the methylation regulation mechanism of chicken in the response of SE inoculation.

Animal and sample collection
Jining Bairi Chicken, a China local chicken breed used in the current study was provided by Shandong Bairi Chicken Breeding Co., Ltd. (Shandong, China). The S. Enteritidis strain (CVCC3377) used in the current study was purchased from the China Veterinary Culture Collection Center (Beijing, China). The animal trial was performed as described previously [91]. In brief, two-day old SE negative Jining Bairi chickens were randomly divided into two groups and raised in two separated isolators with the same condition. Each chicken in the inoculated group was orally inoculated with 0.3 ml 10 9 colony-forming units (cfu)/ml SE inoculant, and each chicken in the control group inoculated with same amount of sterile phosphate buffer saline (PBS). Twelve chickens in each group were euthanized by cervical dislocation at Genomic DNA extraction and DNA library construction Three individual cecum samples from each of inoculated and control groups at 3 dpi were randomly selected for genomic DNA extraction. In total, 6 genomic DNA samples were extracted using TIANamp Genomic DNA Kit (Tiangen, Beijing, China) following the manufacturer's instructions and stored at -80℃ until further use. The concentration and quality of DNA sample was evaluated using DS-11 Spectrophotometer (DeNovix, US) and gel electrophoresis, respectively. Three genomic DNA samples in each of inoculated and control groups were mixed in equal amount to generate a pooled sample for sequencing. The DNA samples were sheared with Covaris ultrasonicator (Life Technology, US). The fragmented DNA was puri ed using AMPure XP beads and end repaired. After end repair and adenylation, Cytosine-methylated barcodes were then ligated to sonicated DNA. Subsequently, 100-300 bp insert size

Bioinformatics analysis of DNA methylation sequencing data
The adaptors and low-quality sequences were removed from the raw data. The clean data were aligned to the chicken reference genome (Gallus_gallus-5.0/galGal5) using the software Bismark according to guideline (coverage > 2X, FDR < 0.05). Uniquely mapped reads were used to estimate the coverage of bases on the whole genome and for further analysis. In addition, we analyzed the genome coverage of the CG, CHG and CHH sites under different sequencing depths, distribution of clean reads in different CG density regions using Bismark and MOABS software [55]. The genome coverage of the CG, CHG and CHH sites in every chromosome and different genome components including promoter, gene body, downstream was also analyzed with the same methods.
According to the localization results of unique mapped reads in the reference genome, the differentially methylated regions (DMR) was detected between the two libraries using Bisul ghter [54] with depth of more than 4X, at least 3 differential methylated sites, the methylation level difference was more than 0.2 (CG type more than 0.3). P < 0.05 was considered as signi cance. The GO enrichment and KEGG pathway analysis was conducted for differentially methylated genes using the BLAST Functional Annotation Tool [92][93][94]. P < 0.05 was considered as signi cance.
Bisul te sequencing PCR analysis (BSP) The speci c primers for BSP were designed using MethPrimer 2.0 (http://www.urogene.org/methprimer2) and listed in Table 4. The bisul te conversion of chicken genome DNA Methylation-Gold Kit following the protocol of EZ DNA Methylation-Gold Kit (Zymo Research, Irvine, CA, USA). The ampli cation converted DNA using Ex Taq Hot Start Version (Takara Bio Inc., Otsu, Japan). And the PCR product was cloned into the pMD18-T vector (Takara Bio Inc., Otsu, Japan). Twenty clones were selected for each fragment and the positive clones were sequenced using ABI3730XL DNA Analyzer (Applied Biosystems, CA, USA). All the sequences were analyzed using BiQ Analyzer [95].

Consent for publication
Not applicable.

Competing Interests
The authors declare that they have no competing interests.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.    The distribution of DMR in different genes. Different color is different location, promoter is the 2000bp of upstream in gene.

Figure 5
The proportion of DMR in each chromosome length.

Figure 6
The distribution of hyper-and hypo-methylated genes in each chromosome   KEGG pathway annotation of differentially methylated genes.

Figure 10
The CpG methylation patterns of HOXD12. Each row is one clone. The black dot is methylated site, the white one is unmethylated site.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.