In our previous studies (Gao et al., 2018a; Gao et al., 2018b), we detected the serum antibody levels for MAP of 8,214 Chinese Holstein cows with the ELISA method, of which 784 cows were identified as positive. Stool samples of each seropositive cow were further collected for quantitative PCR (qPCR) to detect whether MAP was present with an INgene q ParaTB Kit (Ingenasa, Madrid, Spain). Individuals who were positive in both the serum and stool samples and had obvious diarrhea were defined as clinical cows (CC); individuals who were seropositive, but their stools were negative and with no diarrhea, were referred to as subclinical cows (SC); individuals that were negative in both serum and stools samples were treated as healthy cows (HC). Three cows were included in each group (Supplementary Figure S1A). These nine cows were dissected in the same slaughterhouse. The jejunum tissues were collected from each individual within 30 min after slaughtering and placed in liquid nitrogen. The infected part of the jejunum from clinical cows was obviously thickened, inflamed, and edema was present, and that from subclinical cows showed mild inflammation and edema.

Genomic DNA of the jejunum tissue from the nine cows with three in each group was extracted using a TIANamp Genomic DNA Kit (Tiangen, Beijing, China). The DNA quantity and quality were determined using a NanoDrop2000 spectrophotometer (ThermoFisher, Waltham, MA, United States) and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, United States), respectively.

First, 5 μg of lambda DNA was added into the genomic DNA of each sample as the negative control. Then, the genomic DNA was randomly broken into 200–300 bp fragments using an S220 focused-ultrasonicator (Covaris, Santa Clara, MA, United States); these fragments were end-repaired, added to the sequencing adapter, and treated with bisulfite. The methylated cytosines were thereby unchanged and the unmethylated cytosines became uracils, which were changed to thymines after PCR amplification (Supplementary Figure S1B). A library quality control was performed with an Agilent 2100 Bioanalyzer.

With the completion of the library, paired-end sequencing was performed on the Illumina NovaSeq6000 sequencing platform (Illumina, CA, United States). Preliminary quality control of raw reads was carried out with fastqc (v0.11.9, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and these reads were then filtered with fastp software (v0.20) to remove adapters and low quality sequences (Chen et al., 2018). With bowtie2 (v2.4.0), clean reads were deduplicated and aligned against the bovine reference genome (ARS-UCD1.2, https://asia.ensembl.org/index.html), which was bisulfite-converted using bismark (v0.22.1) (Langmead and Salzberg, 2012; Langmead et al., 2018), and the command of “bismark--genome_folder ref.fa -bowtie2 -N 0 -L 20”. Methylation status was determined and methylated CpG sites were marked using the command “bismark_methylation_extractor” (Krueger and Andrews, 2011). The correlation coefficient within each group was calculated using R programing language (v3.6.0).

The bsseq (v1.24.4) package in the R programing language (v3.6.0) was used to find differentially methylated regions (DMRs) among the three comparative groups based on the information about CpG sites (Hansen et al., 2012). Each DMR had at least five methylated CpG sites where the distance between CpG sites was less than 300 bp and there was a greater than 10% difference in methylation levels between groups. The significant ranking of DMRs was dependent on the absolute value of areaStat that was used to determine the significance threshold. In order to screen differentially methylated genes (DMGs), we retained the genes in which the DMRs were entirely located using the “intersect–f 1” command from bedtools (v2.28.0) (Quinlan and Hall, 2010), eliminating genes that had partial DMRs, all these screened genes were defined as DMGs. The promoter is located kilo bases upstream away from the transcriptional start site in the regulatory sequence controlling gene expression (Hernandez-Garcia and Finer, 2014; Yu et al., 2016). Therefore, we defined 2000 bp upstream of the gene as the promoter region and identified the differentially promoter-methylated genes (DPMGs) with identical screening. Both the files of the gene body and promoter were downloaded from UCSC (https://genome.ucsc.edu/). The genes containing both hypermethylated and hypomethylated regions were considered as hypermethylated and hypomethylated genes.

The clusterProfiler tool (v3.16.0) in R was applied to perform GO and KEGG enrichment analysis on DMGs and DPMGs (Yu et al., 2012). The significant threshold related to enrichment analysis of DMGs and DPMGs was a q value <0.05 and p value <0.05, respectively. The p value is a measure of statistical significance; the q value is a measure of the false discovery rate (FDR) (Benjamini and Hochberg, 1995; Storey, 2003).

Genes that were both DMGs and DEGs according to the DNA methylation and RNA sequencing data were identified. Further, we detected the methylated status of promoters in these genes and calculated the association between the methylated status and the fragments per kilobase million (FPKM) with Pearson correlation analysis (Schober et al., 2018).

By performing WGBS, the raw data of the nine jejunum tissue samples from clinical, subclinical, and healthy groups, each comprising three Holstein cows were completed. After quality control, a total of 2.84 × 10⁹ clean reads were obtained with an average of 3.16 × 10⁸ (average 108.33 G per sample) for each sample, indicating a sequencing depth of 30×. There were 81.03% of non-duplicated clean reads uniquely aligned to the bovine reference genome (ARS-UCD1.2, Supplementary Table S1). Bisulfite conversion efficiency (BCE) reached 99.33% suggesting the reliability of the methylome sequencing in this study (Table 1).

We found that on average 3.22% of all genomic C sites were methylated across the nine samples. Methylation in advanced mammals generally exists in three sequence contexts: CpG, CHG (where H is A, C, or T), and CHH. Here, we observed the overall genome-wide methylation levels of 70.50% CpG, 0.40% CHG, and 0.38% CHH in the clinical group, 69.60% CpG, 0.39% CHG, and 0.38% CHH methylation in the subclinical group, and 68.03% CpG, 0.39% CHG, and 0.38% CHH methylation in the healthy group (Supplementary Table S1). Further, we calculated the correlation of 0.9 within each group using the methylation level of CpGs.

To characterize the distribution of methylation in different genomic elements, we analyzed the average DNA methylation levels in the promoter, 5’ untranslated region (UTR), exons, introns, and the 3’ UTR of each methylated gene, and found the average methylation levels of the three experimental groups were similar in the promoter (29.3–30.48%) and 5’ UTR (29.02–30.29%) but displayed certain differences in partial exons, introns, and 3’UTR (Figure 1). Of note, the methylation levels of partial exons, introns, and the 3’UTR in the clinical group (68.50, 72.40, and 73.35%) were obviously higher than those in the healthy groups (67.90, 70.60, and 71.22%). In addition, the methylation level in the X chromosome (63.45%) was significantly lower than that in the autosomes (69.65%, p value <0.05; Figure 2).

By comparing the cytosine methylation profiles of MAP-infected and healthy groups, we detected 100,758 differentially methylated regions (DMRs), including 31,750 between the clinical and subclinical groups, 35,082 between the clinical and healthy groups, and 33,926 between the subclinical and healthy groups (Supplementary Table S2). The length of DMRs was 6–2,292 bp with an average of 350 bp, 98.5% of which were less than 1,000 bp (Supplementary Figure S2). Further, a total of 6,653 differentially methylated genes (DMGs) and 1,779 differentially promoter-methylated genes (DPMGs) were identified, i.e., 3,911 DMGs and 661 DPMGs in clinical vs. subclinical groups, 4,336 DMGs and 893 DPMGs in clinical vs. healthy groups, and 4,094 DMGs and 746 DPMGs in subclinical vs. healthy groups (Supplementary Table S2). Details of the top 10 significant DMGs and DPMGs in the three comparisons are described in Table 2 and Table 3, respectively.

To further investigate the potential associations of the DMGs and DPMGs with paratuberculosis, we implemented functional enrichment analysis with the clusterprofiler package and found that the DMGs were significantly enriched in 77 GO terms (q < 0.05; Figure 3, Supplementary Table S3). From those, the most enriched terms included plasma membrane part (GO:0044459, q = 3.51E-06), phospholipid binding (GO:0005543, q = 3.80E-05), metal ion transmembrane transporter activity (GO:0046873, q = 0.00018), and phosphatidylinositol binding (GO:0035091, q = 0.00223).

According to the KEGG analysis, the DMGs were significantly enriched in 181 pathways (q < 0.05; Supplementary Table S3) with 95 common pathways across all the comparisons (Supplementary Figure S3). These were mainly related to Th1 and Th2 cell differentiation (bta04658, q = 0.03), wnt (bta04310, q = 0.00056), TNF (bta04668, q = 0.0008), MAPK (bta04010, q = 6.98E-08), ECM-receptor interaction (bta04512, q = 0.017), cellular senescence (bta04218, q = 0.03), calcium (bta04020, q = 1.80E-05), and the chemokine signaling pathways (bta04062, q = 0.001).

As for DPMGs, no significant GO terms were enriched (p > 0.05). However, a total of 92 pathways were enriched (p < 0.05, Supplementary Table S3), in which the most significant pathways included oxytocin (bta04921, p = 1.46E-05), influenza A (bta05164, p = 9.65E-05), prostate cancer (bta05215, p = 0.00014), vascular smooth muscle contraction (bta04270, p = 0.00015) and the long-term potentiation signaling pathway (bta04720, p = 0.00022).

Our previous RNA sequencing on the same jejunum samples of the nine Holstein cows used in the present study detected 134 DEGs including 23 between the clinical and subclinical groups, 71 between the subclinical and healthy groups, and 73 between the clinical and healthy groups (q<0.05; unpublished data). Of these, 31 functional genes were differentially expressed and differentially methylated among the three comparative groups simultaneously (Supplementary Table S4). Eight genes, calcitonin receptor like receptor (CALCRL), tenascin C (TNC), GATA binding protein 4 (GATA4), CD44 molecule (CD44), transglutaminase 3 (TGM3), C-X-C motif chemokine ligand 9 (CXCL9), C-X-C motif chemokine ligand 10 (CXCL10), and peroxisome proliferator activated receptor gamma (PPARG), were significantly enriched in immune-related pathways such as Th1 and Th2 cell differentiation, wnt, TNF, MAPK, ECM-receptor interaction cellular senescence, calcium, and chemokine signaling pathways. Therein, CXCL10 and CALCRL contained two methylated sites in the promoter (Supplementary Figure S4). The Pearson correlation coefficients between the methylated status of the promotor and the mRNA levels (FPKM) of CXCL10 and CALCRL were calculated as -0.07 and -0.56, respectively, indicating downregulation roles of the methylated promoter on gene expression (Supplementary Table S5).

Further, nuclear factor of activated T cells 1 (NFATC1), whose methylation of the intron was shown to affect the expression of IL-2 and IL-4 in the effector T cells in previous studies (Akimzhanov et al., 2008; Huang et al., 2011; Christie and Zhu, 2014), was differentially methylated between the MAP-infected and healthy cows and significantly enriched in the Th1 and Th2 cell differentiation pathway (q = 0.033).

Consequently, we inferred the potential of a regulatory mechanism underlying paratuberculosis in dairy cattle. As shown in Figure 4, the differentially methylated NFATC1, TGM3, and GATA4 directly affected the expression of immune factors, including IL-5, IL-2, IL-4, IL-4Ra, IL-2Ra, GATA3, IL-6, IL-12Rb, and IFN-γ secreted in Th1 and Th2 cells, to maintain an inflammatory response. NFATC1 elevated the expression of these immune factors through a cell differentiation pathway. TGM3 activated IL-6 expression via calcium signaling pathways and GATA4 promoted IL-5, IL-4, and IFN-γ through a cellular senescence signaling pathway. PPARG regulated NFATC1 expression through MAPK and calcium signaling pathways thereby modulating immune factors in Th1 and Th2 cells. CXCL9 and CXCL10 stimulated chemokine receptor proteins that promoted the JAK2/3 and STAT pathways, ultimately activating immune factors such as IL-12Rb, IFN-γ, and T-bet. Membrane proteins encoded by TNC and CD44 modulated the expression of IL-5 and IL-4 through the ECM-receptor interaction signaling pathway and another member protein encoded by CALCRL stimulated IL-5 secretion through activation of cAMP in the vascular smooth muscle contraction pathway.

In sum, invasion of MAP into the animal body induced immune responses whereby the above-mentioned nine functional genes stimulated the secretion of immune factors against MAP, indicating the potential regulatory mechanism underlying paratuberculosis in dairy cattle.