Whole genome DNA methylation of triplicate samples of IPEC-J2 cells infected with ETEC F4ab, F4ac, and uninfected, were analyzed to determine methylation profiles of ETEC infection. Using a sliding-window peak-finding algorithm provided by NimbleScan v2.6 (Roche-NimbleGen), a total of 46,940 methylated enrichment peaks (EPs) were identified from the nine samples, of which 14,805 (31.54%) were in the ETEC F4ab infected samples, 16,336 (34.80%) in the ETEC F4ac infected samples, and 15,799 (33.66%) in the uninfected control samples (Table S1). As shown in the methylation map (Figure 1), while most chromosomal regions were covered by methylated peaks, the methylation density in these regions were distinct among the chromosomes; chromosome 13 in particular, contained a relatively large unmethylated region.

In methylome studies, CpG islands (CGIs) are of particular interest because of their role in controlling gene expression (Jones, 2012). Therefore, we analyzed the methylation status of CGIs in the genome of the porcine intestinal epithelial cell line IPEC-J2 after infection with ETEC F4ab and F4ac. We grouped the CGIs into three classes according to their distance to the RefSeq genes: promoter CGIs [from about−10 kb to + 0.5 kb around the transcription start site (TSS)], intragenic CGIs [from + 0.5 kb around the TSS to the transcription terminal site (TTS)], and intergenic CGIs (those that do not fall into neither the promoter northe intragenic group) (Yu et al., 2014; Song et al., 2016; Figure 2A). The numbers of methylated EPs in the three classes of CGIs among the three groups of IPEC-J2 cells are shown in Figure 2B. Most of the methylated EPs were distributed in the intergenic and promoter CGIs for all three IPEC-J2 cell groups. It is worth mentioning that the ETEC F4ab-infected cells had relatively lower methylation levels in promoter CGIs than the uninfected or the ETEC F4ac-infected cells. The ETEC F4ac-infected cells had slightly higher methylation levels in promoter and intergenic CGIs than the uninfected or the ETEC F4ab-infected cells. In addition, the number of methylation EPs within intragenic CGIs in ETEC F4ab/ac-infected cells was higher than in uninfected cells.

Promoters are the main targets of DNA methylation modification, and the center to turn on or off gene expression (Wang X. et al., 2019). The key promoter region is defined as the −800 to +200 bp around the TSS. Based on the CpG ratio, GC content, and length CpG-rich region, we divided gene promoters into three types: high CpG density promoter (HCP), low CpG density promoter (LCP), and intermediate CpG density promoter (ICP) (Yu et al., 2014). We then analyzed the distribution of methylated EPs in the three types of the promoter (Figure 2C). We found that the number of methylated EPs was highest in the HCPs of all IPEC-J2 cell groups, followed by that in ICPs and LCPs. Additionally, in promoter CGIs, HCPs had more methylated EPs in all IPEC-J2 cell groups than ICPs or LCPs (Figure 2D). The activity of HCPs is negatively correlated with their DNA methylation status (Weber et al., 2007). We found that compared to the uninfected cells, the ETEC F4ab infected cells had fewer methylated EPs in HCPs, while the ETEC F4ac infected cells had slightly more methylated EPs in HCPs, this held true for promoter CGIs as well (Figures 2C,D).

To more finely map the pattern of DNA methylation, we divided each type of promoter (HCP, ICP, LCP) into distal (−800 to −200bp) and proximal regions (−200 to +200bp) relative to the TSS. Each region was then defined as methylated (indicated by “1”) or unmethylated (‘0'). Thus, each type of promoter was classified into three subtypes according to their methylation profile (Figure 3A) (Koga et al., 2009). For each type of promoter, the methylation profile and the number of methylated EPs for each IPEC-J2 cell group are graphed in Figures 3B–D. In general, the promoters were either highly distally methylated (“10” pattern) or highly proximally methylated (“01”), but not both (“11”). Infection with ETEC F4ab and F4ac resulted in heterogeneous methylation patterns. For all promoter types, ETEC F4ac infected cells were hypomethylated with ‘01' and “10” patterns compared to uninfected cells. ETEC F4ab infected cells were hypomethylated for the HCP type in all patterns (Figure 3B). For the ICP and LCP types, F4ab infected cells were hypermethylated in the “10” pattern and hypomethylated in the “01” pattern. In the “11” pattern, ICPs, in F4ab infected cells were hypermethylated, and for LCPs they were hypomethylated (Figures 3C,D).

The DNA methylation changes and differences presumably reflect the ability of these pathogens to trigger epigenetic responses in cells involved in the immune system (Tarakhovsky, 2010).

To detect changes in the DNA methylome induced by ETEC F4ab/ac infection, we compared ETEC F4ab/ac infected and uninfected cells to identify differentially methylated regions (DMRs) in genomic DNA. 3297 DMRs were identified between ETEC F4ab infected and uninfected cells, of which 1261 were hypermethylated and 2036 hypomethylated in the infected cells (Table S2). 1593 DMRs were identified between ETEC F4ac infected and uninfected cells, of which 1140 were hypermethylated and 453 were hypomethylated in the infected cells (Table S3). The DNA methylation levels in the DMRs and the differences among the nine samples are shown by using a heatmap in Figure 4A. We subsequently mapped all DMRs to their nearest genomic features and found that DMRs were unevenly distributed across the genome. The majority of the DMRs were in the promoter regions in cells infected with ETEC F4ab or F4ac (Figure 4B), which constitute a large proportion in MeDIP-chip. To predict the potential functional significance of the identified DMRs, we analyzed the DMRs located in intragenic and promoter regions. The DMRs in the F4ab infected cells involve 2601 genes, of which, 2433 were differentially methylated only in promoters, 147 only in gene-bodies, and 21 in both promoters and gene-bodies (Table S4). Of the 2601 genes, 1043 were hypermethylated and 1558 hypomethylated. The DMRs in F4ac infected cells involve 1085 genes, of which, 991 were differentially methylated only in promoters, 91 only in gene-bodies, and 3 in both promoters and gene-bodies (Table S5). Of these 1085 genes, 752 were hypermethylated and 333 hypomethylated. Furthermore, 406 differentially methylated genes were common to both F4ab and F4ac infected cells. Of these, 187 were commonly hypermethylated and 167 commonly hypomethylated, and 52 hyper- or hypo-methylated in opposite directions in F4ab vs. F4ac infected cells (Figure 4C). Some of the identified genes involved in the immune system have been well-studied. For example, PUMA [Bbc3; p53 upregulated modulator of apoptosis (Okuda et al., 2017)] interacts with antiapoptotic Bcl-2 family members, leading to the formation of the free-type Bax and/or Bak, which are then able to signal apoptosis to the mitochondria (Yu and Zhang, 2008). RasGrf1 participates in the Ras signaling pathway (Manyes et al., 2014). JAK2 is in many ways the prototypical member of the JAK family, with an essential signaling role for cytokines and interferons involved in immunity and antiviral responses (Ferrao et al., 2018).

To assess the accuracy of the MeDIP-chip results, the TLR5 gene was selected to validate promoter DNA methylation enrichment using bisulfite sequencing. The TLR5 gene was chosen because results from MeDIP-chip data indicated that it was methylated at low levels in untreated IPEC-J2 cells, moderately methylated in ETEC F4ac infected cells, and highly methylated in ETEC F4ab infected cells (Figure 5A). The results from bisulfite sequencing showed that TLR5 exhibited hypomethylated, moderate methylated, and hypermethylated enrichment in untreated IPEC-J2 cells, ETEC F4ac infected cells, and ETEC F4ab infected cells, respectively (Figure 5B), which are in good agreement with the results from the MeDIP-chip.

We used an Agilent Porcine Oligo Microarray (4 × 44K) to analyze global gene expression in infected and uninfected cells. Genes were classified as differentially expressed if they exhibited |FC (fold-change)| > 1.5 with q < 0.05. In the ETEC F4ab infected cells we identified 606 genes that were significantly differentially expressed compared to the uninfected controls, 421 of these were up-regulated and 185 down-regulated (Table S6). In the ETEC F4ac infected cells, we identified 780 genes that were significantly differentially expressed compared to the uninfected controls, 524 of these were up-regulated and 256 down-regulated (Table S7). This genome-wide expression analysis provided a comprehensive portrait of the immune response, at the transcriptional level, of IPEC-J2 cells infected with ETEC F4ab/ac infection.

To investigate the biological significance of the differentially expressed genes, we performed functional annotation analysis by using a Bioconductor bioinformatics resource. The enriched GO terms could be roughly grouped into two clusters (Figure 6A). The first cluster is factor activity, such as receptor ligand activity, receptor regulator activity, cytokine activity, and chemokine activity. The second cluster centers on receptor binding, e.g., cytokine receptor binding, growth factor receptor binding, chemokine receptor binding, microtubule binding, fatty acid binding, double-stranded RNA binding, and cofactor binding. Since a key step in ETEC infection of pigs is its binding to receptor(s) on enterocytes (Roubos-van et al., 2017), the receptor binding related genes we identified are likely to be involved in this process. KEGG pathway analysis revealed that genes were significantly enriched in some immune/defense response-related pathways (Figure 6B), e.g., IL-17 signaling pathway, TNF signaling pathway and p53 signaling pathway, which indicate their roles in response to the cytotoxic effects of ETEC. Hub nodes have been found to play important roles in networks (He and Zhang, 2006), several hub pathways and associated genes were identified for ETEC F4ab and F4ac infection independently (Figure 7). The genes induced by F4ab and F4ac share most hub pathways, which is consistent with the similarity of the F4ab and F4ac fimbrial antigens. Both have the “a” epitopes formed by the conserved region of the major F4 fimbrial subunit FaeG (Verdonck et al., 2004).

DNA methylation occurring at gene promoters is usually involved in inhibiting the expression of the corresponding genes (Jones, 2012). By using integrative analyses of the DNA methylation profiles and gene expression profiles from IPEC-J2 cells infected with ETEC F4ab/ac vs. uninfected cells, we explored the relationship between methylation changes at the promoter regions and gene expression changes. The results showed that among the differentially methylated genes in cells infected with ETEC F4ab and F4ac, 27 and 10 were inversely expressed with regard to their methylation status, respectively (Tables 1, 2). Genes including S100A9, SGO1, and ESPL1, which are essential for cells' antiviral response, were all hypomethylated and up-regulated in cells infected with ETEC F4ab (Figure 8), these genes were also detected in hub pathways of our transcriptome analysis (Figure 7A). The methylation and expression of MAP3K21, PAK6, and MPZL1, which play roles in immune response and adhesion, were negative correlated in ETEC F4ac infected cells P < 0.05 (Figure 9). Thus, these genes could be powerful candidates methylation target genes related to susceptibility or resistance to ETEC F4ab/ac.

IPEC-J2 cells were grown in Dulbecco's modified Eagle's medium (DMEM)/Ham's F-12 medium (1:1) (GIBCO, Invitrogen, Beijing) supplemented with 5% fetal calf serum (FCS, GIBCO, Carlsbad, CA, USA) and incubated in a humidified 5% CO₂ atmosphere at 37°C. ETEC F4ab strain 195 (O8:K87:F4ab) and ETEC F4ac strain 200 (O149:K91:F4ac) were removed from cryo-storage and cultured in Ordinary Broth Agar at 37°C for three generations (24 h per generation). For cell infection experiments, strains were subcultured in LB medium and incubated with shaking (230 rpm) for 12 h at 37°C. Bacteria were collected by centrifuged and washed with sterile PBS (pH 7.4). Finally, the bacterial suspension was prepared with a final concentration of 1 × 10⁸ CFU/mL (Zhou et al., 2012).

Dose and time of infection optimizations experiments were conducted as described in our previous study (Zhou et al., 2012). Briefly, monolayers of IPEC-J2 cells prepared in 24-well cell culture plates (Corning, Lowell, MA, USA) were washed twice with PBS, then covered with 0.5 mL of DMEM per well. Bacterial suspension was added to the medium of experimental wells (MOI = 10:1) and an equal volume of PBS was aliquoted to the control wells. Each experimental treatment was conducted in triplicate. The cells were incubated for 3 h at 37°C in a 5% CO₂/air atmosphere (Vieira et al., 2010) then collected for nucleic acid isolation.

Genomic DNA and total RNA were exacted by using the Qiagen DNeasy Tissue kit (Qiagen, Hilden, Germany) and TRIzol reagent (Thermo Fisher Scientific, Waltham, U.S.A.) according to the manufacturer's guidelines. Concentrations of DNA and RNA were quantified with Qubit Fluorometer (Thermo Scientific, U.S.A.). Nucleic acid integrity was further analyzed by 1% agarose gel electrophoresis with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA.). For RNA samples, only those with an RNA integrity number >7.0 were retained for subsequent analysis.

Genomic DNA from each sample was sonicated to produce random fragments 200 to 1000 bp in length. Methylated DNA was immunoprecipitated by using BioMag^TM magnetic beads (Bangs Laboratories. Inc) coupled to mouse anti-5-methylcytidine monoclonal antibody (Diagenode). The immunoprecipitated DNA was eluted from the beads and purified by phenol/chloroform extraction and ethanol precipitation. The input DNA and the immunoprecipitated DNA were labeled with Cy3- and Cy5-fluorophere, respectively, and hybridized to an Arraystar Custom Pig CpG Promoter array. The CpG array covered all known CpG islands annotated by UCSC and all Ensembl gene promoter regions (~−800 to +200 bp from transcription start sites), which totally covered ~385,000 probes. Scanning was performed with an Axon GenePix 4000B microarray scanner, following the manufacturer's guidelines detailed in the NimbleGen MeDIP-Chip protocol (NimbleGen Systems Inc., Madison, USA).

Promoter CGIs are defined as CGIs located in-−10 kb to + 0.5 kb around the TSS; intragenic CGIs are located in + 0.5 kb around the TSS to the TTS; the remaining CGIs are defined as intergenic CGIs. Key promoters are defined as the −800 to +200 bp regions around the TSS. The three categories of key promoters were determined as follows: HCP (high CpG density promoter) contains a 500bp interval with CpG/expected CpG above 0.6 and G+C content above 0.55; LCP (low CpG density promoter) do not contain a 500bp interval with CpG/expected CpG above 0.4; and ICP (intermediate CpG density promoter) is the remainders that do not fall into either HCP or LCP (Yu et al., 2014).

The enrichment intensity was determined for each probe that mapped to gene promoters and CGIs in the MeDIP DNA and input DNA samples. Comparisons were represented as log₂ ratio value, where the ratio was defined as the fluorescence signal of MeDIP DNA/the fluorescence signal of input DNA. P-values were calculated to assess whether intensity differences were significant. To avoid technical variability and to evaluate methylation differences between samples, the log₂-transformed ratios were subjected to median centering, quantile normalization, and linear smoothing by using the Bioconductor packages Ringo (Toedling et al., 2007), Limma (Ritchie et al., 2015), and MEDME (Pelizzola et al., 2008). The normalized chip data was then analyzed by using a sliding-window (750 bp) peak-finding algorithm provided by NimbleScan v2.6 (Roche-NimbleGen). A one-sided Kolmogorov-Smirnov (KS) test was applied to determine whether the probes were drawn from a significantly more positive distribution of intensity ratios than those from the rest of the array. Each probe was associated with a –log₁₀[p-value] score from the windowed KS test. If several adjacent probes were significantly above a set threshold, the region was classified as an enrichment peak (EP). NimbleScan detects peaks by searching for at least 2 probes above a –log₁₀[p-value] minimum cutoff of 2. Peaks within 500 bp of each other are merged. The differential probe-level ratio (log₂[MeDIP/Input]) between treated and untreated cells was used to analyze the MeDIP hybridization to identify differentially methylated regions (DMRs) as follows: (1) log₂[ratio treatment]-log₂[ratio control] > 0.25; (2) log₂[ratio control] < 0.4; (3) P-value determined for treatment vs control < 0.05.

Bisulfite sequencing PCR primers was designed by the online MethPrimer software (Li and Dahiya, 2002), which are listed in Table S8. Genomic DNA from the untreated IPEC-J2 cells, ETEC F4ac infected cells and ETEC F4ab infected cells were treated with bisulfite using the EZ DNA Methylation Kit (Zymo Research, D5006). Each experimental treatment was conducted in triplicate. PCR reactions were performed using the ZymoTaq™ PreMix (Zymo Research, E2004) following the manufacturer's protocols. The PCR products were then purified and cloned into the pMD19-T vector (TaKaRa biotechnology Co., Dalian, China). Ten positive clones for each subject were randomly selected for sequencing (Sangon, Shanhai, China). The final sequences were processed by the online software QUMA (Kumaki et al., 2008).

Detailed cDNA microarray procedures are described in our previous work (Zhou et al., 2012). In brief, total RNA was extracted from nine samples (three ETEC F4ab infected cells, three ETEC F4ac infected cells and three control cells) with TRIzol reagent (Thermo Fisher Scientific, Waltham, U.S.A.). Total RNA was amplified and labeled by using a Low Input Quick Amp Labeling Kit, One-Color (Agilent Technologies, Santa Clara, CA, USA), following the manufacturer's instructions. The labeled cRNA was purified by using a RNeasy mini Kit (QIAGEN, GmBH, Germany), then hybridized to porcine oligo microarray slides (Agilent Technologies) containing 43,603 oligonucleotide probes. The hybridized microarray slides were scanned by using an Agilent Microarray Scanner (Agilent Technologies, Santa Clara, CA, USA). Raw data were normalized by using the quantile algorithm from Gene Spring Software 11.0 (Agilent Technologies, Santa Clara, CA, USA). Differentially expressed genes (DEGs) were defined as those with a fold change (|FC|) > 1.5 and P < 0.05.

Gene Ontology (GO) term analyses and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for DEGs were performed by using the R package “clusterProfile” (Yu et al., 2012), with p-values calculated by using right-sided hypergeometric tests. To prevent a high false discovery rate (FDR) in multiple testing, q-values were also estimated for the FDR control (Storey, 2003). Figures were prepared by using the R package “ggplot2” (Wilkinson, 2011). RCircos (Krzywinski et al., 2009) was used to visualize the methylation maps of the nine samples.