A highly virulent strain of serovar 5 (G. parasuis SH0165 strain), was isolated from the lungs of a commercial pig showing signs of arthritis, fibrinous polyserositis, hemorrhagic pneumonia, and meningitis (19). The SH0165 strain was cultured in tryptic soy broth (Difco Laboratories, Franklin Lakes, NJ, USA) or tryptic soy agar (Difco Laboratories) at 37°C, supplemented with 10 µg/mL of nicotinamide adenine dinucleotide (NAD) (Sigma-Aldrich Corporation, St. Louis, MO, USA) and 10% fetal bovine serum (Gibco, Carlsbad, CA, USA).

This study was conducted in strict accordance with the recommendations of the Chinese Regulations for the Administration of Affairs Concerning Experimental Animals 1988 and the Hubei Regulations for the Administration of Affairs Concerning Experimental Animals 2005. The protocol was ratified by the Hubei Province Science and Technology Department of China (permit number SYXK[ER] 2010-0029). Animal studies were approved by the Animal Care and Use Committee of Wuhan Polytechnical University, Hubei Province, China (EM947, 5^th November 2020). All experimental animals were euthanized at the end of the experiment.

Six 28-day-old naturally farrowed early weaned piglets (Duroc × Landrace × Large white), weighing 8-10 kg, were obtained from Wuhan Wannianqing Animal Husbandry Co., Ltd (Wuhan, China) for in vivo experimentation. All experimental piglets were negative for antibodies against G. parasuis when tested using INgezim Haemophilus 11.HPS.K.1 (INgezim, Spain).

The six piglets were randomly divided into two groups including the negative control group and infected group. The infected group was challenged intraperitoneally with 2 × 10⁹ CFU of the G. parasuis SH0165 strain in 1 mL of normal saline. The negative control group was injected intraperitoneally with an equivalent volume of saline alone. All the piglets were monitored for 7 days after challenge, and the cerebrum tissues from the infected and control groups were collected for whole-genome bisulfite sequencing and RNA-seq.

Genomic DNA from porcine brain tissue was extracted with Allprep DNA/RNA/Protein Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. Then extracted DNA (1 µg) was spiked with 1 ng of fragmented unmethylated phage λ DNA (Promega, Madison, WI, USA, average size: 300 bp), which was end-repaired, A tailed, and adapter-ligated using Truseq DNA sample prep kit (Illumina Inc., USA) according to the manufacturer’s protocol. Adapter-ligated DNA was gel isolated and recovered using a QIAquick gel extraction kit (Qiagen, Hilden, Germany). After a cleanup with AMPure XP beads (Beckman Coulter Inc., Fullerton, CA, USA), bisulfite conversion was performed using EZ-DNA Methylation-Gold Kit (Zymo Research, Orange, CA, USA) according to the manufacturer’s protocol. Then, the bisulfite converted DNA was amplified with PfuTurbo Cx Hotstart DNA Polymerase (Agilent Technologies, Inc., Santa Clara, CA, USA) using the following thermal cycles: 95°C 5 min, 98°C 30 s, 12 cycles of (98°C 10 s, 65°C 30 s, 72°C 30 s), 72°C 5 min. Libraries were cleaned with AMpure XP beads and fragment length determined using an Agilent 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA). Quantified libraries (concentration > 2 nM) were then sequenced using the Illumina HiSeq2500 platform.

Total RNA was extracted using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) in accordance with the manufacturer’s instructions and RIN number assess for RNA integrity using an Agilent Bioanalyzer 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA). Quantified total RNA was further purified using an RNAClean XP Kit (Beckman Coulter Inc., Fullerton, CA, USA) and RNase-Free DNase (Qiagen, Hilden, Germany). Only RNA with a RIN of ≥ 7.0 and an 28S/18S ratio of ≥ 0.7 were selected for deep sequencing. Libraries were generated using the VAHTS Total RNA-seq Library PrepKit for Illumina (Vazyme Biotech Co., Ltd., Nanjing, China) and were subsequently sequenced using the Illumina HiSeq X-Ten platform.

General quality control was performed using FastQC v 0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Low quality reads (error ratio > 0.01 or read length < 70 bp or unpaired-reads) and adaptors were removed using Trimgalore v0.4.1 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and Cutadapt v1.9.1 (20). Filtered datasets were aligned to the reference genome, Sus-Scrofa.Sscrofa 10.2, using Bismark v 0.15.0 (21), and using Bowtie2 v 2.2.6 (22) as the underlying alignment tool. Mapping for all datasets generated from the same library were merged, and duplicates removed via the Bismark deduplication tool. Mapped reads were then separated by genome (Sscrofa 10.2 or phage λ) and by source strand (plus or minus). The first four and last one base of each read2 in all read pairs was clipped due to positional methylation bias, and any redundant mapped bases due to overlapping mates from the same read pair were trimmed to avoid bias in the quantification of methylation status. Read pairs mapped to phage λ were used as a QC assessment to confirm that the observed bisulfite conversion rate was >99%. Read pairs mapped to the Sscrofa 10.2 reference genome were used for downstream analysis.

The R package DSS (Dispersion shinkage for sequencing data, DSS) was used to identify DML (Differential methylation loci, DML) with the callDML function, and DMR (Differential methylation region, DMR) with the callDMR function. DMRs that contain at least three validated CpG sites and a minimum length of 50 bp were retained for downstream analysis. DMLs and DMRs with a threshold p value < 0.05 and mean methylation difference > 0.1 were considered to be significantly differentially methylated.

The genomic context of the identified DMLs and DMRs with respect to known genes and gene features was analyzed by using the R package goldmine (http://jeffbhasin.github.io/goldmine). The function getCpgFeatures, getFeatures, and getGenes were used to obtain CpG islands, ENCODE supertracks, and RefSeq gene models from the database of the UCSC Genome Browser (23, 24).

The raw reads were filtered by using the Seqtk sequence processing tool (https://github.com/lh3/seqtk) and mapped to the reference genome (Sus-scrofa.Sscrofa10.2.dna.toplevel.fa) by using the Hisat2 alignment program (25) and StringTie was used to count the gene fragments (26, 27). Package edgeR (28) was used to identify differentially expressed genes with a threshold fold change > 2 and adjusted p-value < 0.05, which was considered as significantly differentially expressed.

To reveal the correlation between DNA methylation status and gene expression, we integrated DNA methylation and gene expression data. First, genes that were significantly differentially methylated in the promoter regions were selected with a threshold of mean methylation difference > 0.1 and p value < 0.05. Second, genes that were significantly differentially expressed were selected with a threshold of |log₂FC| >0.5 and p value < 0.05. For these two selected gene sets, they were firstly merged by the merge function in R. Then, the correlation for the methylation difference and expression difference were determined by “cor. test” in R using Pearson method. Significantly negative correlations were considered if the correlation coefficient was negative and the p value was < 0.05, and significantly positive correlations were considered if the correlation coefficient was positive, and the p value was < 0.05. Differentially expressed genes that were negatively correlated with methylation were subjected to further studies.

To discover the potential biological function underlying the identified genes, we conducted GO and KEGG enrichment analysis on the differentially expressed or differentially methylated genes. Genes were mapped to each term within the GO database. GO terms with adjusted p values < 0.05 were considered as significantly enriched. Pathway annotation of the entire genome was conducted by the KEGG automatic annotation server. Pathways with p values ≤ 0.05 were deemed to be significantly enriched.

STRING database (ELIXIR, Hinxton, Cambridge shire, UK) were utilized for the protein-protein interaction (PPI) network construction, based on the methylation and expression of inversely correlated genes. Network Analyzer (http://www.mpi-inf.mpg.de/) was used to determine the topological properties of the biological networks. Cytoscape software (Institute for Systems Biology, Seattle, WA, USA) was used for network visualization.

Bisulfite sequencing was performed to validate the DMRs detected by the whole-genome DNA methylation analysis. DNA (330 ng) was isolated using Quick-DNA Miniprep Kit (Zymo Research, Orange, CA, USA) and DNA bisulfite conversion was conducted by using the EZ DNA Methylation Gold Kit (Zymo Research, Orange, CA, USA). Bisulfite-converted DNA (1 µg) was amplified by PCR with EpiTaq polymerase (for bisulfite-treated DNA) (TaKaRa Biotechnology Co., Ltd., Dalian, China), and primers ( Supplementary File 1 ) were manually designed using Primer Premier 5. Before primer design, DNA sequences need to be methylation converted in the MethPrimer website (http://www.urogene.org/methprimer/). PCR products were cloned using the Topo TA cloning kit (Invitrogen Corporation, Carlsbad, CA, USA) and after verification of successful cloning, at least ten different clones for each sample were sent for sequencing at TSINGKE Biological Technology. The DNA methylation profile of individual CpG sites were analyzed and presented using BiQ Analyzer (https://biq-analyzer.bioinf.mpi-inf.mpg.de/).

Total RNA was extracted from porcine brain tissues using RNeasy Mini Kit (Qiagen, Hilden, Germany), and reverse transcribed into cDNA using reverse transcriptase (TaKaRa Biotechnology Co., Ltd., Dalian, China). cDNA was further quantified using the SYBR Green PCR Kit (TaKaRa Biotechnology Co., Ltd., Dalian, China) according to the manufacturer’s protocol. Individual samples underwent three technical repeats, with β-actin as the reference gene. The Primer Premier 5.0 software was utilized to design the PCR primers ( Supplementary File 2 ).

Experimental data were analyzed using the two-tailed Student’s t-test and presented as the mean ± standard deviation and differences with P < 0.05 were considered to be significant (^* P < 0.05; ^** P < 0.01).

Whole-genome bisulfite sequencing was performed with the Illumina HiSeq 2500 to explore the global DNA methylation changes associated with the pathology of meningitis induced by G. parasuis. A total of 626.50 ± 47.40 million and 627.76 ± 36.65 million raw reads were obtained from the infected and negative control groups, respectively. After filtration and quality control, the clean reads ratio ranged from 96.1% to 99.0% ( Supplementary File 3 ). All the multi-mapped reads duplicated reads and adaptors were removed and aligned to the Susscrofa genome, with the unique mapping ratio ranging from 54.3% to 60.7% ( Supplementary File 3 ). Unmethylated phage λ DNA was added to evaluate the bisulfite conversion ratio, which ranged from 99.1% to 99.4% ( Supplementary File 3 ).

Cytosine methylation levels at a single base-pair resolution were explored using callDML function in the R package DSS. Results showed that 316,621 CpG sites were differentially methylated in G. parasuis infected porcine brain, with 132,557 CpG sites hypermethylated and 184,064 CpG sites hypomethylated (P < 0.05) ( Supplementary File 4 ). We also identified 12,145 (DMRs) that were statistically significant between the two groups, with 4,649 DMRs hypermethylated and 7,496 DMRs hypomethylated in G. parasuis infected porcine brain (P < 0.05) ( Supplementary File 5 ).

The genomic and CpG island context of the DMRs were investigated in G. parasuis infected porcine brain and results showed that the differentially methylated regions were distributed throughout the whole genome ( Figure 1 ). DMRs were most common in intergenic regions (64.7%), and less common in introns (26.1%). Only a minor percentage of DMRs occurred in promoters (3.9%), followed by exons (2.8%), transcriptional termination sites (TTS, 1.6%), 3’ untranslated regions (3’ UTR, 0.6%) and 5’ untranslated region (5’ UTR, 0.4%) ( Figure 1A ). The genomic annotation profile exhibited similarities when we focused on hyper- and hypo- methylated DMRs separately. The majority of DMRs were annotated in intergenic regions (64.4% and 64.9% for hypermethylated and hypomethylated DMRs, respectively), and a minor percentage of DMRs were annotated in promoters (4.7% and 3.3% respectively) ( Figure 1C ). Furthermore, the CGI (CpG island) annotation showed that most DMRs were found outside CGIs and most of these were in the open sea regions (greater than 4 kb from CpG islands) ( Figure 1B ). For hypermethylated and hypomethylated DMRs, 14.7% and 10.0% were annotated in CGIs, 12% and 11.2% were in shores (2000 bp flanking CpG islands), and 5.4% and 5.7% were in shelves (2000 bp flanking shore regions) ( Figure 1D ).

RNA-seq analysis was performed to explore changes in the global transcriptome occurring in meningitis induced by G. parasuis infection. Unsupervised hierarchical clustering identified two unique clusters that had distinct signatures ( Figure 2A ). A total of 1,126 genes were significantly differentially expressed in G. parasuis infected porcine brain, including 606 up-regulated and 520 down-regulated genes (fold change > 2, adjusted p < 0.05) ( Figure 2B and Supplementary File 6 ).

GO and KEGG enrichment were performed to reveal the potential function of the affected genes involved in the pathogenesis of meningitis induced by G. parasuis. GO enrichment showed that cellular processes, single-organism processes, and metabolic processes were the most abundant terms in the biological process category. In the cellular component category, cell, cell part and membrane were the top three terms, while binding, catalytic activity, and molecular transducer activity were the top terms in the molecular function category ( Figure 2C ). KEGG analysis revealed that cytokine-cytokine receptor interaction, cell adhesion molecules (CAMs), and phagosome were the most significantly abundant pathways involved in the pathogenesis of G. parasuis induced meningitis ( Figure 2D ).

DNA methylation and gene transcription are normally causal related, where hypermethylation might be associated with gene down-regulation and hypomethylation might be related to gene up-regulation, especially when methylated sites are located in the promoter region of genes. The hierarchical clustering of the top-10000 significantly differentially methylation sites, considering different genomic context (5’UTR, 3’UTR, exon, intergenic, intron, promoter-TSS, and TSS), showed a striking difference between the control and G. parasuis infected samples ( Figure 3A and Supplementary Figure 1 ). A volcano plot showed the distribution of different methylation sites under different methylation thresholds ( Supplementary Figure 2 ). To further identify the significance of DNA methylation in the regulation of gene expression, we integrated the CpGs methylation patterns with the gene expression profiles obtained from the same sample. The differentially methylated CpG sites in the promoter region and the differentially expressed mRNA transcripts were aligned to chromosomes throughout the whole genome ( Figure 4 ). The integrated analysis showed that 301 differentially expressed genes corresponding to 1068 differentially methylated CpGs in their promoter regions exhibited an inverse correlation (|log₂FC| > 0.5, |diffMethy| > 0.1, and P < 0.05) ( Figure 3B and Supplementary File 7 ). Furthermore, 141 genes (28%) were significantly hypermethylated and down-regulated, while 160 genes (29%) were significantly hypomethylated and up-regulated (|log₂FC| > 0.5, |diffMethy| > 0.1, and P < 0.05) ( Figure 3B ). The inversely correlated methylation sites were detected throughout the genome on all chromosomes except chromosome 19, 20 and chromosome Y ( Figure 3C ).

The same integrated analysis was performed on the DMRs and gene expression profiles. Forty DEGs corresponding to 40 DMRs in the promoter regions were inversely associated, and 33 DEGs were positively correlated with their promoter DMRs (|log₂FC| > 0.5, |diffMethy| > 0.1, and P < 0.05) ( Figure 3D and Supplementary Files 8, 9 ). Table 1 shows the top 30 significant DEGs with an inverse correlation between DMR and gene expression.

In addition to the inversely correlated genes, we also observed an equivalent number of positively correlated genes, which were hypermethylated but up-regulated or hypomethylated but down-regulated ( Figures 3B, D ). Thus, both negative and positive correlations between DNA methylation and mRNA expression were observed in our G. parasuis infected porcine brain samples, suggesting that there may be two different mechanisms involving DNA methylation-dependent gene regulation for G. parasuis induced meningitis.

GO and KEGG enrichment were performed to reveal the biological dysfunctions caused by the DNA methylation and the expressed inversely correlated genes in G. parasuis infected samples. GO enrichment showed that cellular processes, single-organism processes, and metabolic processes were the most abundant terms in the biological process category. In the cellular component category, cell, cell part, and organelle were the top three terms, while binding, and catalytic activity were the top terms in the molecular function category ( Figure 5A ). When compared to the KEGG pathways enriched by all the differential genes ( Figure 2D ), cell adhesion CAMs, and hematopoietic cell lineages were common between the two enrichments, and bacterial invasion of epithelial cells, and RIG-1-like receptor signaling pathways were particularly enriched by the inversely correlated genes in G. parasuis infected porcine brain ( Figure 5B ).

To determine the reliability of the sequencing data, mRNA expression and DNA methylation at base-pair resolutions were examined on two inversely correlated genes, semaphoring 4D (SEMA4D) and von Willebrand factor A domain containing 1 (VWA1), according to the significance criteria. Bisulfite sequencing showed that the methylation rate in the DMR of SEMA4D decreased from 65.4% to 39.1% after G. parasuis stimulation, while the methylation rate in the DMR of VWA1 decreased from 45.8% to 21.5% after G. parasuis stimulation ( Figures 6A, B ). Validation with RT-PCR showed that the mRNA expression of SEMA4D and VWA1 were significantly up-regulated after G. parasuis stimulation (P < 0.01), consistent with the hypomethylation in the promoter regions of each gene ( Figure 6C ).

To reveal the potential associations within the differentially methylated genes, we created a protein-protein interaction network which consisted of 193 protein nodes and 243 interactions based on the differential genes (|log₂FC| > 0.5, |diffMethy| > 0.1, and P < 0.05) ( Figure 7A ), with CCL5 and CXCL10 as the core genes. The crucial roles of genes were determined by the node degree and the top 40 genes with the highest connection degrees were selected and subsequently used for KEGG pathway enrichment ( Figure 7B ). The KEGG result showed that the top 40 genes mainly participated in chemokine signaling pathways, cytosolic DNA-sensing pathways, RIG-1-like receptor signaling pathway, Toll-like receptor signaling pathway, CAMs signaling pathways, and TNF signaling pathway ( Figure 7C ).