Roving methyltransferases generate a mosaic epigenetic landscape and influence evolution in Bacteroides fragilis group

Abstract

Three types of DNA methyl modifications have been detected in bacterial genomes, and mechanistic studies have demonstrated roles for DNA methylation in physiological functions ranging from phage defense to transcriptional control of virulence and host-pathogen interactions. Despite the ubiquity of methyltransferases and the immense variety of possible methylation patterns, epigenomic diversity remains unexplored for most bacterial species. Members of the Bacteroides fragilis group (BFG) reside in the human gastrointestinal tract as key players in symbiotic communities but also can establish anaerobic infections that are increasingly multi-drug resistant. In this work, we utilize long-read sequencing technologies to perform pangenomic (n = 383) and panepigenomic (n = 268) analysis of clinical BFG isolates cultured from infections seen at the NIH Clinical Center over four decades. Our analysis reveals that single BFG species harbor hundreds of DNA methylation motifs, with most individual motif combinations occurring uniquely in single isolates, implying immense unsampled methylation diversity within BFG epigenomes. Mining of BFG genomes identified more than 6000 methyltransferase genes, approximately 1000 of which were associated with intact prophages. Network analysis revealed substantial gene flow among disparate phage genomes, implying a role for genetic exchange between BFG phages as one of the ultimate sources driving BFG epigenome diversity.

Fig. 1. Four-decade collection of clinical BFG isolates from the NIH Clinical Center.

MLST marker gene cladogram of BFG genomes sequenced in this study supplemented with Genbank reference genomes (n = 462 total). Taxonomy assignments were defined proteomically with MALDI-TOF mass spectrometry (Bruker Biotyper) and genomically with GTDB-Tk. “Source” and “Decade” data were extracted from clinical laboratory metadata records.

Fig. 2. Genomic analysis demonstrates open pangenomes and substantial movement of mobile genetic elements within and between species.

a Stacked barplots quantifying average numbers of persistent, shell, and cloud genes across eight species. b Pangenome analysis for a subset of BFG species. Rarefaction curves indicate open pangenomes over the sequenced set, with the three largest pangenomes demonstrating greater than 20,000 genes each. c Analysis of accessory region and mobile genetic element content. Top panel shows species-level spread bins for more than 33,000 accessory regions/mobile genetic elements. “Species” indicates number of species sharing the indicated number of accessory regions or mobile genetic elements. Paired bottom panel barchart indicates annotated features of the accessory regions as a percentage of each paired spread-level bin in the upper panel.

Fig. 3. DNA methyltransferases are remarkably diverse and abundant in the accessory genomes of BFG species.

Host species are represented across the top row as squares (area proportional to number of genomes analyzed). DNA methyltransferase gene families (80% AAI, 80% AF) are represented in rectangular grids below as filled circles (with area proportional to number of genes in the family) with colors indicated in the key. Edges connect species with DNA methyltransferase gene families that are encoded by one or more genomes within the species. Location of the given methyltransferase gene family in the core, shell or cloud genome in is indicated by edge color, and edge thickness indicates the number of times the gene family is encoded in the genome of the species. In this analysis, ‘Core’ was defined as presence in >90% of genomes in a species, ‘Shell’ was defined as presence in >10% and ≤90% genomes in a species, and ‘Cloud’ was defined as presence in <10% of genomes in a species.

Fig. 4. Network analysis reveals substantial methyltransferase gene flow among disparate phage genomes.

Network graph of phage viral operational taxonomic units (vOTU) clusters and DNA methyltransferase gene families (80% AAI, 80% AF), indicated by nodes of different shapes as defined in the legend. Shape size is proportional to the number of phage genomes within a given vOTU cluster or methyltransferase genes within a gene family. Edges connect methyltransferase gene families and vOTU clusters containing prophage genomes that encode a methyltransferase gene from that methyltransferase gene family. Edge thickness is proportional to the number of genomes that encode the corresponding gene family.

Fig. 5. BFG methylation motifs tile a vast combinatoric space with open panepigenomes over the sequenced isolate set.

a Rarefaction curves of DNA methylation motifs detected in genomes of BFG species in this study and comparison with C. difficile and B. breve species from external studies (denoted with ‘*‘; data^,). BFG rarefaction curves indicate open panepigenomes over the sequenced isolate set. b–f Heatmaps of detected DNA methylation motifs in Parabacteroides distasonis, B. vulgatus, B. fragilis (sensu stricto), B. ovatus, and B. thetaiotaomicron isolates. Rows indicate individual isolates with corresponding marginal MLST marker gene phylogenies. Columns indicate distinct methylation motifs. The three most prevalent motifs in each set are labeled, and labels are omitted for the rest. Cells are colored when a given motif is present in the corresponding isolate and colors indicate the class of base modification as indicated in the legend. Sequences of all motifs can be found in Supplementary Data 10.

Fig. 6. DNA methylation motifs demonstrate genome-wide depletion in lineage-specific context.

(Top) A swarm plot of Z scores of density (motifs/kilobase) of each motif across genomes. Two-sided T-tests were performed for each target motif along with its control motifs and Benjamini–Hochberg testing performed separately with FDR < 1%. Unadjusted p values are reported as follows: *p < = 0.05, **p < = 0.01, ***p < =1e-3, ****p < =1e-4; “ns” indicates that the reported p value was determined to be non-significant after testing with Benjamini-Hochberg at FDR = 1%. (Middle) Data displayed as motifs per kilobase. Boxes show data quartiles and whiskers show 1.5X IQR with diamonds representing outliers. (Bottom) MLST marker gene cladograms by species using all methylome-analyzed genomes with heatmap of lineage-specific DNA methylation motifs. (All) B. fragilis genomes, n = 108. B. ovatus genomes, n = 44. B. vulgatus genomes, n = 28.