Hemin availability induces coordinated DNA methylation and gene expression changes in Porphyromonas gingivalis

ABSTRACT Periodontal disease is a chronic inflammatory disease in which the oral pathogen Porphyromonas gingivalis plays an important role. Porphyromonas gingivalis expresses virulence determinants in response to higher hemin concentrations, but the underlying regulatory processes remain unclear. Bacterial DNA methylation has the potential to fulfil this mechanistic role. We characterized the methylome of P. gingivalis, and compared its variation to transcriptome changes in response to hemin availability. Porphyromonas gingivalis W50 was grown in chemostat continuous culture with excess or limited hemin, prior to whole-methylome and transcriptome profiling using Nanopore and Illumina RNA-Seq. DNA methylation was quantified for Dam/Dcm motifs and all-context N6-methyladenine (6mA) and 5-methylcytosine (5mC). Of all 1,992 genes analyzed, 161 and 268 were respectively over- and under-expressed with excess hemin. Notably, we detected differential DNA methylation signatures for the Dam “GATC” motif and both all-context 6mA and 5mC in response to hemin availability. Joint analyses identified a subset of coordinated changes in gene expression, 6mA, and 5mC methylation that target genes involved in lactate utilization and ABC transporters. The results identify altered methylation and expression responses to hemin availability in P. gingivalis, with insights into mechanisms regulating its virulence in periodontal disease. IMPORTANCE DNA methylation has important roles in bacteria, including in the regulation of transcription. Porphyromonas gingivalis, an oral pathogen in periodontitis, exhibits well-established gene expression changes in response to hemin availability. However, the regulatory processes underlying these effects remain unknown. We profiled the novel P. gingivalis epigenome, and assessed epigenetic and transcriptome variation under limited and excess hemin conditions. As expected, multiple gene expression changes were detected in response to limited and excess hemin that reflect health and disease, respectively. Notably, we also detected differential DNA methylation signatures for the Dam “GATC” motif and both all-context 6mA and 5mC in response to hemin. Joint analyses identified coordinated changes in gene expression, 6mA, and 5mC methylation that target genes involved in lactate utilization and ABC transporters. The results identify novel regulatory processes underlying the mechanism of hemin regulated gene expression in P. gingivalis, with phenotypic impacts on its virulence in periodontal disease.

T eeth are secured in the mouth through the periodontal ligament, which is attached to the subgingival root surface and the alveolar bone. Periodontal disease is a common chronic inflammatory disease of the periodontal tissues that is driven by bacteria in subgingival microbial biofilms on the root surface. In periodontal disease, the integrity of the attachment tissues becomes compromised leading to bleeding, gingival recession, formation of deep periodontal pockets, and periodontal bone loss, ultimately culminating in tooth loss in a susceptible host. Periodontitis is a global health issue affecting both developed and developing nations (1). Furthermore, periodontitis is also associated with a range of extra-oral inflammatory conditions, such as diabetes and cardiovascular disease, rheumatoid arthritis, Alzheimer's disease, and pregnancy complications (2)(3)(4).
During the progression from health to periodontal disease, the oral microbiota undergoes a major transition wherein the microbial community increases in total bacterial diversity, accompanied by a bloom in disease-associated bacteria that start to dominate the community, while being otherwise present in low numbers in healthy microbiota (5,6). Porphyromonas gingivalis, a black-pigmenting anaerobic rod bacterium, is well established as one of the organisms, which is positively associated with this shift to a dysbiotic microbiota in disease. Indeed, P. gingivalis has been shown to act as a trigger for the conversion of the subgingival microbiota to a more diverse and more pathogenic community in animal models of disease (5,6). This keystone pathogen behavior has been attributed to the expression of a variety of extracellular virulence determinants including those trafficked through a type IX secretion system (T9SS) and through the production of large numbers of outer membrane vesicles, which are collectively able to destabilize the host-microbe homeostasis at this mucosal surface (7,8).
Hemin, iron (III) protoporphyrin IX chloride, is a major source of iron for a variety of organisms, including P. gingivalis, where it is an essential co-factor for growth. Porphyro monas gingivalis has a dedicated transport system to ensure an efficient accumulation of hemin on the surface for iron utilization. Subsequent conversion of hemin to μ-oxo bishaem, and its complexation with iron III, leads to black pigmentation of the bacterial cells, a prominent characteristic of P. gingivalis on solid media containing blood (9). Porphyromonas gingivalis HmuY and HusA are dedicated, surface-exposed, hemophorelike proteins that scavenge and bind hemin, free or released from heme-containing proteins by gingipain proteases. The complexes deliver the bound heme to outer-mem brane receptors (HmuR or HusB) for internalization, and heme is then imported through a dedicated TonB-dependent transport machinery (10).
The concentration of environmental hemin acts as a regulator of virulence in this bacterium: cells cultured in high concentrations of hemin-comparable to the situation in a diseased periodontium with elevated bleeding-upregulate the expression of a wide variety of virulence determinants, including a range of hydrolytic enzymes such as the gingipain cysteine proteases, and other protein products, many of unknown function (11,12). High-hemin grown cells are significantly more pathogenic in animal models of tissue destruction compared to cells grown under hemin limitation (7,(13)(14)(15). Hence, the availability of hemin in the growth media is critical to the full expression of virulence of P. gingivalis, and hemin may be regarded as a global regulator of transcription in this bacterium. Although previous studies have characterized the impact of hemin availabil ity on gene expression and virulence in P. gingivalis, the regulatory processes underlying these functional changes are not well understood. However, the P. gingivalis ATCC33277 genome encodes seven two-component signal transduction system (TCS) pairs of a response regulator and a histidine kinase traditionally involved in gene regulation in other organisms (16). Recently, it has been reported that one of these TCS, PorX/PorY, responds to environmental hemin and, through the co-ordinated action of PorX and SigP, regulates transcription of the T9SS and hemin accumulation; mutants defective in either porX or porY are attenuated in a mouse model of virulence (17,18). However, given the magnitude of changes in response to the concentration of hemin, there may be other mechanisms of hemin-regulated gene expression in P. gingivalis.
Like eukaryotic genomes, bacterial genomes are subject to epigenetic modifica tions, specifically DNA methylation, which has multiple roles in bacteria including the regulation of transcription. Bacterial DNA methylation includes modifications N6-methyladenine (6mA), 5-methylcytosine (5mC), and N4-methylcytosine (4mC), of which 6mA is most prevalent in prokaryotes. Recent development in single-molecule sequencing technologies, for example, Nanopore sequencing, allows for the detection of these multiple DNA methylation base modifications. Using this approach, bacterial DNA methylation signatures have recently been shown to impact gene expression profiles and genome stability (19,20).
In this study, we investigated the impact of hemin availability on the P. gingivalis epigenome and transcriptome. Porphyromonas gingivalis W50 was grown in controlled conditions, in continuous culture at a constant pH and identical growth rates, under either limited hemin or excess hemin conditions previously associated with virulence changes in the same system (21). Subsequently, DNA sequencing using Nanopore and RNA sequencing using Illumina were carried out to profile the P. gingivalis DNA methylome and transcriptome under each condition. Hemin-dependent differentially modified signatures were identified for gene expression, as well as for 6mA and 5mC DNA methylation profiles both in Dam/Dcm motifs and genome wide. Comparison of the signals highlighted a cluster of coordinated changes in six genes, including in genes related to lactate utilization and in ABC transporters.

Continuous culture chemostat configuration
A 1-L, sealed glass chemostat was autoclave sterilized with 300 mL BHI in situ. One hundred milliliters of P. gingivalis culture was inoculated into the chemostat through a designated inoculation tube, using a peristaltic pump (Watson Marlow 101U). Sterile BHI (supplemented with either 0.2 or 5 mg L −1 hemin for limited and excess conditions, respectively) was continuously fed at a flow rate of 20 mL h −1 to achieve the working volume of 700 mL, before being increased to 70 mL h −1 , in conjunction with har vest/waste pump activation at the equivalent rate. These conditions set a fixed dilution rate of 0.1 h −1 , which corresponds to constant mean generation time of 6.9 h. Steady state conditions were achieved after approximately 10 mean generations. Environmen tal conditions were controlled/monitored using an Applikon ADI 1030 Bio-controller and BioXpert V2 software (Applikon Biotechnology, The Netherlands). Environmental conditions were controlled for included nutrients, temperature, oxygen, and pH. The pH of the culture was maintained automatically at pH7±0.1 using 1 M NaOH and 0.5 M HCl fed through a peristaltic pump. Temperature was preserved at 37 ± 0.3°C with a silicon heat jacket. The chemostat was continuously sparged with filtersterilized 5% CO 2 /95% N 2 to maintain anaerobic conditions and stirred at 40 rpm (Motor Controller ADI 1012, Applikon Biotechnology; Supplemental Note, Figure A).

Experimental design
The dose-response curve of P. gingivalis' exposure to hemin was previously determined by us in similar chemostat conditions (15). McKee et al. (15) showed that hemin amounts greater than 1 mg L −1 increased concentration of hemin in the chemostat without impacting cell yield. Because growth medium is added continuously at a fixed rate in the chemostat, hemin will always be either "limiting" growth or present in "excess. " Here, two independent chemostat experiments were performed, one with 5 mg L −1 hemin (excess) and one with 0.2 mg L −1 hemin (limitation).
Post-inoculation, the optical density (OD 600 ) of the culture was monitored every 12 h using a benchtop spectrophotometer (Jenway 6305). Steady state was defined as three consecutive OD 600 readings within 10% of each other. Total viable counts and culture purity were determined daily, following serial dilutions in phosphate buffered saline, inoculation of CBA plates in triplicate, and growth in a Whitley A45 anaerobic Worksta tion (Don Whitley Scientific, UK). Sampling for DNA and RNA analysis was performed every 48 h during steady state (Fig. 1). This represented biological replicates separated by five full-volume replenishments of the chemostat vessel. Overall, three biological replicates were collected for each condition, with each one of the replicates collected at one sampling point (days 10, 12, and 14).

Nucleic acid sampling
During nucleic acid sampling periods, the harvest bottle was replaced with a sterile collection bottle (on ice) and approximately 20 mL was collected. For DNA sampling, 1 mL aliquots were centrifuged at 13,000 rpm for 10 min (4°C), prior to supernatant removal and storage at −80°C. For RNA samples, ten 1-mL samples were similarly centrifuged, the supernatant decanted, and the pellet reconstituted in 1 mL of DNA/RNA Shield (Zymo Research, USA) prior to storage at −80°C.

DNA extraction
One-milliliter cultures were centrifuged at 17,000×g, at 4°C for 10 min. The pellets were used to extract high molecular weight DNA using Gentra Puregene yeast/bacteria DNA isolation kit (Qiagen, Germany), following the manufacturer's protocol. DNA was finally resuspended in 200-μL EB buffer (Qiagen) and stored at 4°C before use.

Illumina RNA sequencing
RNA samples were first subjected to rRNA depletion using the NEBNext rRNA Depletion Kit (Bacteria) (New England Biolabs, USA). The samples were then used for RNA-Seq library preparation using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina. Both kits were used according to the manufacturer's instructions. Sequencing was performed by Novagene, UK on an Illumina Novaseq6000 instrument.
Oxford Nanopore DNA sequencing MinIon library preparation and DNA sequencing was performed at the University of Nottingham on MinIon flow cell (Oxford Nanopore, UK). Whole-genome sequencing libraries were prepared from six bacterial genomic DNA samples. DNA samples were quantified using the Qubit 4 Fluorometer (Thermo Fisher Scientific, USA) and the Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific; Q32853), and the Agilent 4200 TapeSta tion and the Agilent Genomic DNA ScreenTape Assay (Agilent, USA; 5067-5365 and 5067-5366) were used to assess the molecular weight of the DNA. Sequencing libra ries were prepared using the Native barcoding genomic DNA protocol (Oxford Nano pore Technologies; Version: NBE_9065_v109_revS_14Aug2019). This protocol uses the Ligation Sequencing Kit (Oxford Nanopore Technologies; SQK-LSK109) and the Native Barcoding Expansion 1-12 kit (Oxford Nanopore Technologies; EXP-NBD104) to barcode the unsheared DNA fragments within each sample. Barcoded samples were then pooled in equimolar amounts for the final sequencing adapter ligation. All purification steps were performed using AMPure XP beads (Beckman Coulter; A63882). The final barco ded library pool was loaded onto a MinION flow cell (Oxford Nanopore Technologies; FLO-MIN106 R9.4.1) and run on the GridION X5 Mk1.

RNA sequencing data processing
Sequencing data were processed using the ProkSeq pipeline v2.8 (22) Briefly, FastQC (23) and AfterQC (24) were implemented to check the quality of data, trim adaptor sequen ces, and filter out low-quality reads using standard parameters. This resulted in a total of ≈25.3-36.0 M PE reads (2 × 76 bp) being kept for downstream processing. Subse quently, bowtie2 (25) was used to align reads against the 104 contigs of the P. gingivalis W50 genome assembly deposited in the NCBI RefSeq (accession no. GCF_000271945.1). RefSeq gene annotations were extracted from P. gingivalis' GTF file, and the number of total reads per gene were calculated using featureCounts (26). Counts per gene were estimated for 1,992 genes. RSeQC (27) was used to study sample coverage uniformity along the body of genes.

Motif-independent DNA methylation calling
Nanopore signals were originally stored in multi-read fast5 files, and the ont_fast5_api "multi_to_single_fast5" (https://github.com/nanoporetech/ont_fast5_api) was used to generate single fast5 files before detection of non-standard nucleotides with Tombo v1.5.1 (30). First, electric current signal level data were re-squiggled using the same P. gingivalis W50 reference genome used to preprocess the RNA-Seq and methylated motif data. Secondly, all-context motif-independent models were used in Tombo to detect the alternative bases 6mA and 5mC from the signal data. Alternative basecalling in each sample was summarized per genomic position with coverage and dampened fractions estimated at each adenine and cytosine site of both positive and negative DNA strands. Coverage was stored in the bedGraph format and proportion of methylation in the WIG format. Data files were imported into R using the package rtracklayer (31) and combined into a single data frame using GRanges, IRanges, and findOverlaps from the IRanges package (32).
The proportion of methylation-dampened fraction-was estimated by adjusting the proportion of methylated reads according to coverage at low-coverage sites in order to include them in downstream analyses (30). Despite this, only adenines and cytosines with at least 10× coverage across samples were included in our analysis here as well. Overall, 1,150,659 adenines and 1,079,809 cytosines were kept for downstream analyses, with 350-525× mean base coverage per sample.

Statistical analyses
Differential gene expression analysis DESeq2 (33) was used to normalize the read count data and perform differential gene expression analysis. Growth in limited hemin was the "control" and growth in excess hemin was the "treatment. " All three samples from each condition were used, and global sample differences were observed using the principal component analysis (PCA) method. Both DESeq2 and edgeR (34) were implemented from ProkSeq's downstream analysis. Both methods allow for inter-sample comparison in differential gene expression analysis. Here, we used edgeR results as sensitivity for DESeq2 results since edgeR normalizes expression data differently from DESeq2 (DESeq2 uses the median of ratios method, and edgeR uses the trimmed mean of M values to normalize read count data). The false discovery rate (FDR) method was used to identify differentially expressed genes (DEGs) and genes with FDR <5%, and a log 2 fold change (LFC) >1.5 or <−1.5 were respec tively considered over and under-expressed in excess hemin. DEGs were subsequently manually annotated to the P. gingivalis W83 genome (GenBank accession no. AE015924) (35) in order to compare results with previous research. Synteny analysis at specific loci was performed using MAUVE (36).

Gene ontology enrichment analysis
Gene ontology (GO) enrichment analysis in over and under-expressed gene sets was performed using ShinyGO (37). Because the genome of P. gingivalis W50 was not available from ShinyGO's species database, only well-characterized genes that matched gene symbols in P. gingivalis W83 reference strain genome were used in this analysis. A more relaxed FDR threshold of 10% was applied to assess enrichment in this case and GO molecular functions were reported.

Differential DNA methylation analysis
Mean and median DNA methylation across samples were calculated using all motifs, adenines, and cytosines kept for downstream analysis. Global DNA methylation differences between the samples were observed using PCA for Dam/Dcm and all-con text 6mA and 5mC modifications, selecting for minimum 10× and 100× coverage per position.
Two-sided Student's t-tests were applied at each motif, adenine, and cytosine site after filtering for positions where at least a 5% mean methylation difference was observed between the experimental conditions [79 "GATC" motifs (Dam), 139,794 all-context adenines, and 87,398 all-context cytosines]. A sensitivity analysis was performed for minimum 100× coverage [66 "GATC" motifs (Dam), 85,090 all-context adenines, and 46,930 all-context cytosines] to ascertain whether the significance of results was impacted by depth of sequencing. Results were considered statistically significant at an FDR threshold of 5% to match gene expression analysis.
The differentially methylated motifs and all-context positions identified were annotated to nearby genes considering the location and gene strand. First, differen tially methylated "GATC" motifs ("GATC"-DMMs), differentially methylated adenine sites (DMAs), and differentially methylated cytosines (DMCs) were annotated to the gene they were located if present in the gene body, and secondly, they were annotated to genes found in the same contig and within ±1 kb of their location. Genes located in the opposite strand of the methylation signal were considered in both cases.
In comparing the overlap between DEGs, DMA, DMCs, and "GATC"-DMMs, we considered all DMAs, DMCs, and "GATC"-DMMs that were in or near the gene (±1 kb away, both strands). DMA and DMC data were further used to test the wider DNA methylation effects on the expression of nearby genes. A contingency table contain ing all 1,992 genes was constructed, and Fischer's exact tests were employed to test for the enrichment of genes with concurrent gene expression and all-context DNA methylation changes. The enrichment of differential DNA methylation directly upstream, downstream, and in the body of genes was also tested. Here, Wilcoxon signed-rank tests were used to test for generalized shifts in DNA methylation signals between the experimental conditions at sites 100 bp upstream, 100 bp downstream, and within the body of genes previously annotated to DMAs and/or DMCs.

Motif enrichment analysis
The MEME (Multiple EM for Motif Elicitation) online tool (38) was used to detect enriched motifs containing the DMAs and DMCs identified with concurrent gene expression changes. For this, 15 base pair long DNA sequences surrounding the DMAs and DMCs identified were extracted. Subsequently, motifs with minimum 4 bp and maximum 15 bp were searched using MEME's Classic objective function and the ZOOPS distribution model. Both the given strand and the reverse complement strand were considered, and the top 10 motifs found were reported. Motifs were considered significant if their reported E-value threshold was ≤0.05.

RESULTS
Gene expression and DNA methylation profiling of P. gingivalis W50 was performed after culturing in chemostat runs with excess (5 mg L −1 ) or limited hemin (0.2 mg L −1 ). The dose-response curve of P. gingivalis' exposure to hemin was previously determined by us (15), where we showed that 0.2 mg L −1 of hemin limits P. gingivalis' growth and 5 mg L −1 of hemin results in excess hemin. We chose to use chemostats for these studies because this approach provides cells at a physiological steady state under constant environmen tal conditions. Growth of P. gingivalis in both high and low hemin was at a constant specific growth rate, and all culture parameters remained constant with the exception of the concentration of hemin. Biological replicates were collected every 48 h after achieving steady state, with three independent biological replicates per experimental condition profiled using Illumina RNA-Seq and Nanopore DNA sequencing technologies (Fig. 1). culture [mean = 1.89 × 10 9 colony forming units (CFU mL −1 ), range 5.07 × 10 8 -3.27 × 10 9 CFU mL −1 ] versus limitation (mean = 1.80 × 10 8 CFU mL −1 , range 1.63 × 10 8 -2.05 × 10 8 CFU mL −1 ). OD 600 ranged between 1.67-1.72 (mean 1.69) and 0.357-0.376 (mean 0.367) for excess and limited hemin, respectively. The pH was automatically controlled at pH 7.0 ± 0.1, while temperature remained between 36.7 and 37.1°C for both conditions. Bacteria were grown at an identical and constant growth rate of 6.9 h mean generation time under both conditions (15).

Gene expression changes
RNA-Seq data for the six samples was processed using ProkSeq (22) and DESeq2 (33) to assess gene expression changes observed after culturing P. gingivalis W50 in excess (treatment) and limited (control) hemin conditions. These results were mapped to the P. gingivalis W50 genome (RefSeq assembly accession no. GCF_000271945.1). The expression data captured all 1,992 previously annotated genes, and those with an expression LFC >1.5 or LFC <−1.5 were considered differentially expressed after multiple testing correction at FDR = 5%. Altogether, 161 and 268 DEGs were respectively overand under-expressed in excess hemin across the P. gingivalis genome (Fig. 3).
Different types of genes involved in virulence were identified in both the under-and over-expressed DEG subsets. For example, these included different types of N-acetylmur amoyl-L-alanine amidase genes, T9SS type A sorting domain-containing protein genes and MATE family efflux transporter genes, and tetratricopeptide repeat protein genes (Tables S1 and S2). Additionally, methyltransferase genes were also identified in both the over-and under-expressed DEG subsets (e.g., tRNA and 16S rRNA methyltransferases, and class I SAM-dependent methyltransferases, respectively).
We explored molecular functions shared across the DEGs, although a proportion of genes were not well characterized. Sixty-seven under-expressed (25% of total) and 19 over-expressed (12%) genes identified during growth in excess hemin encoded hypothetical proteins. Additionally, only 16 (6%) of under-expressed genes and 49 (30%) of over-expressed genes identified had well-characterized annotations in RefSeq. A GO term enrichment analysis was carried out using ShinyGO and gene symbols that matched genes in the reference strain P. gingivalis W83 (35), with altogether 8 (out of 16 well-characterized) under-expressed and 40 (out of 49 well-characterized) overexpressed genes. Forty-two and 43 GO molecular functions were respectively enriched in genes under-and over-expressed in excess hemin conditions. In both cases, molecular functions identified were mostly related to binding activities such as cyclic compound binding, heterocyclic compound binding, and ion binding ( Fig. 4; Tables S3 and S4). In excess hemin, functions related to Fe-S cluster binding, 4Fe-4S cluster binding, and metal cluster binding, and several different catalytic activities, including GTPase, Research Article mSystems metallo-endopeptidase, active transmembrane transporter, and diphosphotransferase activities, were also identified ( Fig. 4; Table S4). Two sensitivity analyses were carried out for the genome-wide DESeq2 DEG results presented here. A first sensitivity analysis was performed using edgeR, which normalizes expression data differently from DESeq2, and consistently identified the majority of DEGs identified with DESeq2 (Supplemental Note, Figure B). Specifically, edgeR identified 166 genes as over-expressed and 266 genes as under-expressed with excess hemin, of which 161 and 265 were previously identified with DESeq2. In a second general sensitivity analysis, a PCA of the genome-wide RNA-seq levels also confirmed overall global differences in gene expression levels between the culturing conditions (Supple mental Note, Figure C).

DNA methylation changes
Nanopore sequencing data were generated for the same six biological samples analyzed with Illumina RNA sequencing. Dam-and Dcm-dependent methylation was character ized using the Guppy basecaller. A total of 12,066 "GATC" (Dam) and 819 "CCWGG" (Dcm) motifs with minimum 10× coverage across samples were considered. Motif-independent methylation levels for 6mA and 5mC were characterized using Tombo with "all-context"  Table S1. LFC, log2 fold change; lfcSE, standard error of log2 fold change.
Research Article mSystems models, considering >2M positions with a minimum 10× coverage. PCA of the genomewide DNA methylation levels pointed to both global 6mA and 5mC differences in P. gingivalis according to hemin culture conditions [Supplemental Note, Figure  . Differential methylation analyses with respect to hemin levels were then carried out, considering signals where at least 5% mean methylation difference between the experimental conditions was observed after multiple testing correction (FDR 5%).

Dam/Dcm DNA methylation
DNA methylation signals were analyzed at 79 "GATC" motifs in the genome of P. gingivalis W50, which had at least 5% mean methylation differences between hemin conditions. Of these, 36 motifs showed statistically significant differential methylation signals across hemin conditions [Table S5 (A)]. "GATC"-DMMs were located in alto gether 32 genes, including a zinc-dependent metalloprotease (HMPREF1322_RS05790), a thioredoxindisulfide reductase (HMPREF1322_RS01475), and a 4-alpha-glucanotrans ferase (HMPREF1322_RS03650). If genes up to ±1 kb away from the "GATC" motifs were considered, 85 genes were candidates for possible Dam-like regulation. Sensitivity No minimum 5% mean methylation difference was observed at each Dcm-associ ated "CCWGG" motif, in line with previous observation for no overall differences in the Dcm-based PCA analysis [ Supplemental Note, Figure E(A/B), bottom panels]. This suggests either the absence of enzymes targeting "CCWGG" motifs in P. gingivalis W50 or a minimal role of a Dcm-like response to hemin availability.

All-context 6mA/5mC DNA methylation
The proportion of differential methylation at each genomic position was analyzed for altogether 139,794 adenines and 87,398 cytosines with at least 5% mean difference in DNA methylation levels between conditions. Forty-nine DMAs were identified between experimental conditions across the P. gingivalis W50 genome (Fig. 5). A cluster of 15 DMAs was observed in a genomic region (NZ_AJZS01000011.1 ≈12-20 kb) containing the (Fe-S)-binding protein HMPREF1322_RS00720, lactate utilization protein HMPREF1322_RS00725, hypothetical protein HMPREF1322_RS00730, PaaI family thioesterase HMPREF1322_RS00735, and ABC transporter ATP-binding protein HMPREF1322_RS00740 ( Fig. 5 top panel, purple; Fig. 6). In annotating DMAs to genes, if genes up to ±1 kb away from the DMA on both strands were considered, all DMAs are annotated to at least one gene for a total of 75 genes altogether. Furthermore, 25 DMAs were annotated to genes on the same strand, including seven DMAs encoding a lactate utilization protein (HMPREF1322_RS00725), five in a hypothetical protein (HMPREF1322_RS00730), and two each in a 4-alpha-glucanotransfer ase (HMPREF1322_RS03650) and a Ppx/GppA family phosphatase (HMPREF1322_RS06180) ( Table 3). Four DMAs were located within 1 kb of a gene encoding a (Fe-S)-binding protein (HMPREF1322_RS00720). In a sensitivity analysis considering only adenines with at least 100× coverage (85,090 adenines), 20 DMAs were identified at FDR = 5% (Table S6), and all 20 were also identified in the results of the main analysis.
Forty-seven DMCs were identified between experimental conditions across the P. gingivalis W50 genome (Fig. 5). Two clusters of 15 and 12 DMCs emerged-the first was the same as the DMA cluster region (NZ_AJZS01000011.1, ≈12-20 kb), and the second was in a genomic region (NZ_AJZS01000050.1, ≈0-3.5 kb) containing the DUF2436 domain-con taining protein HMPREF1322_RS04325 and nucleoside permease HMPREF1322_RS04330 ( Fig. 5 bottom panel, purple and green; Fig. 6). Thirtyfive DMCs were located in 16 genes Darker nodes are more significantly enriched; bigger nodes represent larger gene sets; and thicker edges correspond to more overlaps between gene sets. Pathways were created based on gene matches to P. gingivalis W83 reference genome.

Research Article mSystems
including in genes identified in the DMA analyses above. Specifically, four DMCs were located in genes encoding the (Fe-S)-binding protein HMPREF1322_RS00720, two in the lactate utilization protein HMPREF1322_RS00725, and two in the hypothetical protein HMPREF1322_RS00730 reported above for 6mA. Single DMCs were also identified in the 4-alpha-glucanotransferase HMPREF1322_RS03650 and Ppx/GppA family phosphatase HMPREF1322_RS06180 genes (Table 4). Overall, if genes up to ±1 kb away on both strands were considered, 31 genes harbored DMCs, of which 10 also harbored DMAs.
Other genes annotated to the DMCs identified included those encoding transporter proteins such as an iron ABC transporter permease (HMPREF1322_RS01240) and others (HMPREF1322_RS01245; HMPREF1322_RS06460). In sensitivity analysis considering only cytosines with at least 100× coverage (46,930 cytosines), 47 DMCs were also identified, of which 27 were also in the main analysis results (Table S7). We explored evidence for putative DNA sequence motifs at the identified all-context DMAs and DMCs. First, we observed that the DMAs identified in this all-context analysis ( Table 3) did not overlap with our "GATC"-DMM results (Table S5). We observed, however, that the 4-alpha-glucanotransferase HMPREF1322_RS03650 and the DUF3298 domaincontaining protein HMPREF1322_RS07770 contained respectively two and one DMAs within 500 bp of a "GATC"-DMM. Overall, DMAs and "GATC"-DMMs shared 13 genes annotated within 1 kb of their location (out of 75 and 85 genes respectively identified for DMAs and "GATC"-DMMs) (Table S5). Second, we searched for over-representation of up to 15-nucleotide long sequences surrounding all DMAs and DMCs. Although the analyses identified putative sequences of length 10-11 nucleotides, the results did not reach statistical significance (MEME E-value >0.05; Supplemental Note, Figure F).

Genes with coordinated differential DNA methylation and expression changes
We assessed whether the observed DEGs also exhibited DNA methylation changes for 6mA and 5mC. We first compared DEGs to DMAs and DMCs each (Fig. 7), and then explored further genomic regions that contained DEGs, DMAs, and DMCs (Fig. 7).
Overall, 23 DEGs also harbored DMAs (Table 5 and 6). Of these, eight genes were over-expressed in excess hemin including the SDR family oxidoreductase HMPREF1322_RS07305 (P. gingivalis W83: PG2069) and the Ppx/GppA family phospha tase HMPREF1322_RS06180 (P. gingivalis W83: PG1739) genes (Table 5). In all the eight cases, the DMA was on the same strand as the gene, located either in or upstream of the gene (Table 5), with hypermethylation in excess hemin for most   Table 6). Seven of the 15 genes had DMAs located upstream of the gene, and six of these showed increased 6mA accompanied by decreased gene expression, indicating a potential regulatory effect of 6mA in dampening the expression of P. gingivalis genes according to hemin availability in the bacterium's microenvironment.
Most DMAs annotated to under-expressed genes had higher adenine methylation in excess hemin (84%) regardless of their location with respect to gene structure. In comparison to all-context DMAs, the "GATC"-DMMs were annotated to 14 DEGs (six and eight DEGs respectively over and under-expressed; Tables S8 and S9). DEGs associated with "GATC"-DMMs included among others the 4-alpha-glucanotransferase HMPREF1322_RS03650 (W83: PG0767) mentioned immediately above. However, only four genes had "GATC"-DMMs upstream their start site, and another two genes showed greater methylation associated with lower expression overall. This suggests very limited potential for a role of GATC methylation in the regulation of gene expression in P. gingivalis W50 in response to hemin availability. One over-expressed and 11 under-expressed genes also had DMCs within 1 kb of the gene (Table 7 and 8). Two of the 11 under-expressed genes had DMCs upstream of the gene (Table 8) with increased methylation levels observed for all 11 genes. Moreover, 95% of DMCs annotated to under-expressed genes had increased cytosine methylation in excess hemin regardless of their location in relation to the gene. This is consistent with results from 6mA, although knowledge of gene regulation by 5mC methylation in bacteria remains limited.
For all DEGs that also harbored DMAs or DMCs, we analyzed the overall levels of adenine (Supplemental Note, Figure G) and cytosine (Supplemental Note, Figure  H) methylation upstream (−100 bp), downstream (+100 bp), and in the gene body. Differential adenine and cytosine methylation were observed across all regions analyzed for the lactate utilization protein HMPREF1322_RS00725 (W83: PG1172) (upstream; gene body; downstream; Table 9; Fig. 6). Differential adenine and cytosine methylation were also observed in the gene body of the ABC transporters HMPREF1322_RS00740 (W83: PG1175) and HMPREF1322_RS00745 (W83: PG1176) and hypothetical protein HMPREF1322_RS00730 (W83: PG1173) (Supplemental Note, Figures G and H). These ABC transporter proteins are located directly downstream of a lactate utilization gene cluster, which includes the lactate utilization protein and hypothetical protein, possibly highlighting a key region of the P. gingivalis genome regulated by DNA methylation (Fig. 6). Overall changes in adenine and/or cytosine methylation were also, for example, observed in the body of the 4−alpha−glucanotransferase HMPREF1322_RS03650 (W83: PG0767) gene, and upstream the ABC transporter HMPREF1322_RS00740 (W83: PG1175) (Supplemental Note, Figures G and H).

DISCUSSION
We identified a significant impact of hemin availability in the growth medium on the P. gingivalis epigenome and transcriptome. Growth of P. gingivalis is restricted in the healthy gingival crevice due to a tight restriction by the host on hemin availability. However, during increasing inflammation in periodontal disease, the availability of heme-containing molecules is substantially increased. Hemin-dependent differentially modified signatures were identified for gene expression, as well as for 6mA and 5mC methylation profiles. Comparison of the signals highlighted a cluster of coordinated changes in genes related to lactate utilization and in ABC transporters. The findings identify putative bacterial DNA methylation signatures that may play a regulatory role in the expression of genes, in response to hemin availability.
Many of the genes that were over-expressed in excess hemin encoded proteins that contain iron as part of their tertiary structures or have iron/hemin binding domains. These include Lys-gingipain (Kgp; Table S1) that is known to be involved in iron acquisition and virulence (14). Kgp is a large multi-domain protein, which is able to hydrolyze hemoglobin, bind hemin via its hemagglutinin domains, and facilitate its acquisition and conversion into the black pigment (μ-oxo bishaem) at the cell surface (40). Previous work confirmed virulence changes in response to hemin availability under similar experimental conditions (21).
Over-expressed members of the ABC transport complex (Fig. 6B and Table 6), forming a classical operon, may potentially represent a novel hemin or alternative solute acquisition system that may influence the cellular phenotype. The enhanced expression of thiamine phosphate synthase is not only limited to this gene (Table S1) of the pathway; the contig also encodes ThiH, ThiG, ThiC, ThiS, and a hypothetical protein, all arranged in an apparent operon. Thus, the biosynthesis of the energy co-factor thiamine, vitamin B1, may constitute an important aspect of P. gingivalis metabolism. Rubrerythrin offers protection of the anaerobe, P. gingivalis, against the toxicity of molecular oxygen and hydrogen peroxide; an isogenic mutant is labile under these conditions (41). In addition, there appears to be a preference for succinate/fumarate and a-ketoglutarate utilization (Table S1).
Many genes under-expressed in excess hemin were involved in hemin binding or transport. Some genes also encoded for DNA transfer functions such as TraC, TraE, TraF, and TraJ. Interestingly, the gene encoding the largest protein component of the T9SS, Sov (SprA), and (42) also genes encoding OmpH, PorN (GldN), and PorX along with five cargo proteins, were under-expressed in excess hemin (Table S2). Sov and PorN are both major structural proteins of the T9SS located at the outer membrane. PorX is a member of a two-component regulator system (PorX/PorY) that governs transcription of numerous genes including those encoding the T9SS (42). OmpH is a cargo protein of the T9SS able to cause activation of proinflammatory cytokines (interleukin-6 and tumor necrosis factor α) in macrophages (43). The T9SS complex transports over 30 cargo proteins, including known virulence factors, with a CTD signal motif, to the external environment (44). However, the entire structure and organization of this system is not fully elucidated, and it is not currently clear why only a subset of this multi-protein complex and a limited number of the numerous cargo proteins are under-expressed in excess hemin. T9SS has many components with differences in relative stoichiometry. Furthermore, different genomic regions and operons contribute to express components of the T9SS (45), introducing a further level of complexity such as variable response of individual promoters to different environmental exposures.
The gene encoding the long fimbriae surface appendage, FimA, together with its accessory proteins of the fim operon were also under-expressed in excess hemin. The genes encoding regulatory proteins [e.g., AraC, TetR (AcrR), and LuxR] followed a similar expression pattern. Many transposases were also under-expressed in excess hemin, whereas no transposases were over-expressed in excess hemin, suggesting that the expression of these mobile elements is a characteristic of a stress environment where access to iron is limited. Time of incubation in the chemostat can lead to amplified expression changes in some of the genes (e.g., the ABC transporter genes in the lactate utilization cluster; Supplemental Note, Figure I).
Epigenetic mechanisms, including the methylation of specific DNA sequences by DNA methyltransferases, allow unicellular organisms to respond rapidly to environmental stresses (46). Here, we showed that both adenines and cytosines can be differentially methylated in P. gingivalis according to hemin availability in its microenvironment. Unlike    in eukaryotes, in bacteria, 5mC is not the dominant DNA modification and instead exists alongside 4mC and 6mA, the latter being the most prevalent modification in prokaryotes (47). Here, we characterized both 6mA and 5mC, and verified higher levels of adenine methylation (mean = 23% and median = 13% for 6mA genome-wide) in comparison to cytosine methylation (mean = 18% and median = 13% for 5mC genome-wide). However, lack of previous DNA methylation studies in P. gingivalis and a general lack of knowledge of the most common genetic motifs targeted by DNA methyltransferases (DNMTs) in this bacterium limited the analyses that we could undertake, and may impact some of the interpretation of data. We analyzed DNA methylation changes in Dam/Dcm motifs (mean = 8% and median = 13% for Dam; mean = 0.02% and median = 0% for Dcm genomewide) to investigate possible Dam-and Dcm-like DNA methylation in P. gingivalis. However, we observed only limited evidence for Dam-like "GATC" methylation effects that overlapped differential gene expression. Alternatively, and because no standard model for methylation calling in P. gingivalis exists, we used Nanopore reads and allcontext predictive models in Tombo to estimate the proportion of reads harboring methylated bases in our study. All-context models in turn introduced a higher degree of uncertainty in methylation calling in our analysis. Additionally, we used proportion of methylated reads as an estimate of methylation levels in each sample of bacterial cells. In our data, most positions at which DNA methylation was detected were not fully methyla ted or demethylated, pointing to presence of heterogeneity in the level of DNA methyla tion across different bacterial cells from the same culture sample (Supplemental Note, Figures K and L). Heterogenous DNA methylation in subpopulations from a single colony can be responsible for the phenotypic variation observed during a bacterium's adapta tion to a new environment (46). In this case, the DNA methylation levels observed could indicate a flexible response to hemin that is implicated in the colonization of the gingival crevice by P. gingivalis in vivo. However, further work is needed to test the role of methylation heterogeneity in P. gingivalis. DNA methylation has historically been associated with bacterial DNA restriction-mod ification systems, which protect bacterial cells from foreign DNAs (46). More recently, DNA methylation has also been shown to play important roles in other aspects of bacterial biology such as timing of DNA replication, timing of transposition, and conjugal transfer of plasmids, DNA repair, and cell partitioning, and regulation of gene expression (46). Less is known, however, about how widespread the role of DNA methylation is in virulence and pathogenesis. Prior evidence suggests that the activity of DNMTs such as Dam can, for example, affect virulence in a number of pathogens including E. coli, Salmonella spp., Yersinia pseudotuberculosis, and Vibrio cholerae (48).
Here, we showed that the availability of hemin-a critical micronutrient that modulates the expression of virulence factors in P. gingivalis-is linked to changes in the DNA methylome of this bacterium. We also showed that methylation changes were coordinated to alterations in the expression of lactate utilization and ABC transporter genes; however, their exact role in P. gingivalis remains unknown. Porphyromonas gingivalis is assaccharolytic and utilizes amino acids/peptides as the major source of nutrition (49). Hence, any changes in lactate utilization are unconnected to the fermentation of monosaccharides via glycolysis and reduction of pyruvate. Lactate in this bacterium may arise through an endogenous metabolic mechanism via degrada tion of amino acids or via a lactate import system, which has also been shown to be upregulated in hemin excess. Other stressors including oxygen and growth in defined medium can also affect lactate metabolism in P. gingivalis. In Lewis et al. (50), hmuY and a lactate permease (PG1340) were both downregulated in micro-aerophilic conditions; and extracellular lactate was significantly reduced in anaerobiosis, while intracellular lactate was unaffected. In Moradali and Davey (51), exogenously added lactate to bovine serum albumin-based media downregulated virulence-associated genes including sprA (sov), ragA, and fimA, and significantly reduced bacterial growth. Thus, lactate flux is an important aspect of P. gingivalis physiology under different environmental conditions. It should be noted, however, that in the current study, the only variable between     DE analysis, results from the differential expression analysis; Distance Meth, distance of the methylated position to the gene start site in base pairs, negative values indicate upstream location; DM analysis, results from the differential methylation analysis; ExH, growth in excess hemin; LFC, log2 fold change; LiH, growth in limited hemin; Strand Expr, strand of the gene; Strand Meth, strand of the methylated position.
the two conditions was the hemin concentration: all other variables were constant in these continuous culture experiments. Previous research found that lactate utilization (in Haemophilus influenzae and Staphylococcus aureus) and the ABC transport of iron, manganese, and/or zinc ions across the membrane (e.g., in Yersinia pestis, Salmonella enterica and Streptococcus spp.) are linked to increased virulence (52)(53)(54). In our study, differential methylation in lactate utilization and ABC transporter genes co-occurred with their decreased expression in excess hemin conditions. We hypothesize that the decreased expression of these genes in excess hemin, and conversely their relative increased expression in limited hemin, could be related to P. gingivalis scavenging of nutrients according to their availability in the growth environment. Functional assess ment of a methylation-mediated response to hemin is however necessary and will be the scope of future experiments. Porphyromonas gingivalis possesses the AI-2 [encoded by luxS (HMPREF1322_RS03160, PG0498) (55-58); quorum sensing system]. The expression of LuxS is cell density-dependent. James et al. (57) suggested that luxS (PG0498 in P. gingivalis W83) is required for hemin uptake since a luxS isogenic mutant is affected in genes involved in hemin/iron acquisition (57). Furthermore, separate microarray and transcriptome analysis of the luxS mutant (59, 60) indicated the involvement of quorum sensing in hemin utilization in P. gingivalis, specifically through upregulation of hmuY and hmuR in DluxS. In our results, both hmuY and hmuR are under-expressed, to different extents, in excess hemin. It is therefore possible that increased cell density afforded by increased environmental hemin may play a role in the regulation of gene expression observed in this study. Limiting hemin conditions caused a decrease in cell density prior to achieving steady state in our results (Fig. 2). We hypothesize that this could be due to hemin limitation not being optimal for P. gingivalis growth. Alternatively, the starting OD may be skewed by the presence of dead cells in the inoculating culture, where these are later removed in the chemostat run.
In conclusion, using Nanopore sequencing, we were able to characterize for the first time the epigenetic landscape of P. gingivalis W50, and to compare epigenetic variation to changes in gene expression. We showed that iron availability in the bacterium's microenvironment can lead to discrete DNA methylation changes potentially impacting the expression of genes related to growth and virulence in humans. Future work is needed to experimentally dissect the role of DNA methylation on the observed gene expression changes in P. gingivalis. Furthermore, there is a need to better characterize the diversity of motifs targeted by different DNMTs, including from P. gingivalis. A more comprehensive characterization of DNMTs and their targeted motifs in the bacterial world will allow better future exploring of DNA methylation in P. gingivalis with Oxford Nanopore. Methodology, Resources, Supervision, Visualization, Writing -original draft, Writingreview and editing

DATA AVAILABILITY STATEMENT
The Nanopore and Illumina RNA sequencing data generated for this study have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB60468.