Abstract

Main

Fn, a member of the oral microbiota, has gained attention as an emerging cancer-associated bacterium. Worldwide, unbiased genomic analyses have revealed an enrichment of Fn in human CRC relative to non-cancerous colorectal tissues⁹. Previous work by our group and others demonstrated that patients with CRC tumours harbouring high levels of Fn have poorer survival⁵, that Fusobacterium colonizes regions of patient tumours with immune and epithelial cell functions supportive of cancer progression¹⁰, that Fusobacterium persists in metastatic disease⁶ and that microbiome modulation targeting Fn could change the course of this disease^6,11. Moreover, exogenous Fn infection in animal and cellular models has supported a cancer-promoting role for this bacterium^6,9–14. However, considerable strain-to-strain variation in Fn genotypic and phenotypic features has been described^12–16. Such heterogeneity has raised challenges with reproducing Fn cancer-inducing phenotypes in animal and cellular models^16,17 and it has been proposed that only a select group of Fn strains may possess carcinogenic capabilities¹⁷.

Here, leveraging a comprehensive collection of human CRC Fn strains and carrying out extensive comparative genomics, we reveal that a select clade within Fn subspecies predominates the CRC niche. In vitro and in vivo functional studies demonstrate that this clade is highly virulent in the context of CRC, and should be a primary focus in subsequent mechanistic studies on Fn pathogenicity in CRC and for the development of targeted inhibitors.

Niche-enriched Fn genes and subspecies

Classical microbiology culture approaches have re-emerged as valuable tools to functionally assess members of tissue-associated microbiomes. Here we carried out Fusobacterium targeted culture on 130 human CRC tumours from which 65 Fusobacterium CRC-associated strains were obtained from 59 unique patients. Given that Fn is predominantly an oral pathobiont, we included 81 Fusobacterium strains isolated from the oral cavity of individuals without cancer, as a control group. These oral strains were obtained from the American Type Culture Collection (ATCC) and the Korean Collection for Oral Microbiology (KCOM) repositories. Using PacBio long-read single-molecule real-time¹⁸ sequencing, we generated complete and closed genomes, with corresponding epigenetic methylomes, for 146 unique Fusobacterium strains (Fig. ​(Fig.1a1a and Supplementary Tables 1 and 2), 92% of which belonged to the species Fn (n=135; n=55 CRC associated and n=80 oral associated; Fig. ​Fig.1b).1b). As Fn is the most frequently detected species in CRC tumours^2,3, we therefore focused our analysis on a comparison of these 55 CRC-associated and 80 oral-associated Fn genomes.

Open in a separate window

Fig. 1

Fn niche features.

a, A schematic of Fusobacterium strain collection (n=146) and the sequencing strategy for unique strains. SMRT, single-molecule real-time sequencing. b, A column graph depicting the proportion of Fusobacterium genomes, subset by species, within the CRC (orange) and oral (blue) niches. The inset shows all non-Fn species of Fusobacterium (Fnec, F.necrophorum; Fu, F.ulcerans; Fp, F.pseudoperiodonticum; Fc, F.canifelinum; Fv, F.varium). c, The composition of the Fn pangenome subset by niche. Anvi’o²¹ gene cluster (GC) prevalence was used to define core (≥95%), accessory (≥5% and <95%) and rare (<5%) features conserved in both CRC-associated and oral-associated strains (collection core, ≥95% in all strains within the collection; collection cloud, ≥5% and <95% in all strains within the collection; collection rare, <5% in all strains within the collection). Disparate features are those that do not fall into any of the other noted bins. d, The proportion of niche-enriched gene clusters across CRC-associated and oral-associated Fn genomes. The plot box shows the 25th percentile, median and 75th percentile. The plot whiskers indicate the minima and maxima. e, KofamKOALA KEGG orthologue analysis²⁷ of niche-enriched gene clusters. f, A column graph depicting the proportion of Fn genomes, grouped by subspecies, within the CRC and oral niche. Statistical analysis was carried out using a two-sample z-test, two-tailed. NS, not significant. g, Gene presence–absence heat map of canonical Fn virulence factors (fadA (refs.^38,39,42), fap2 (ref.³⁴), fplA (ref.³³), radD (refs.^36,62), aim1 (ref.³⁵), cmpA (ref.³⁷) and fusolisin³²) across Fn subspecies, in which each column represents an individual genome (Fna n=75, Fnn n=17, Fnp n=33, Fnv n=10). The heat map is organized using an rpoB gene-based phylogenetic tree. For each genome, the tree end points indicate the niche origin (CRC (orange); oral (blue)) and the bar colour indicates the Fn subspecies (Fna (red); Fnn (gold); Fnp (purple); Fnv (brown)). The graphics in a were created using BioRender.com.

Given that Fn strains in human CRC tumours are predicted to originate from the human oral cavity^19,20, but are rare members of the lower gastrointestinal (GI) tract microbiota of individuals without cancer¹, we reasoned that CRC-associated Fn strains harbour an additional genetic repertoire to facilitate their colonization in human CRC tumours. To test this, we examined our 135 Fn genomes using the analysis and visualization platform for ‘omics data (Anvi’o) workflow for microbial pangenomics²¹. Pangenomic analysis identifies all genes present in a species (‘pangenome’) and discerns between gene content conserved among most members (‘core genome’) and gene content shared among subsets of members (‘accessory genomes’)^22,23. We observed that accessory genome size increases as the number of sampled Fn genomes increases, supportive of previous proposals that Fn has an open pangenome²⁴ (Extended Data Fig. ​Fig.1a).1a). To account for uneven sampling²⁵ between CRC (n=55) and oral (n=80) genomes, we subset our analysis on the basis of niche and found that CRC-associated Fn strains have a smaller accessory genome compared to that of oral-associated strains (Fig. ​(Fig.1c,1c, Extended Data Fig. ​Fig.1b1b and Supplementary Table 1). Functional enrichment analysis²⁶ identified 483 and 241 gene clusters significantly enriched (q<0.05) in CRC and oral strains, respectively (Fig. ​(Fig.1d1d and Supplementary Table 3). Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologue analysis²⁷ of the 724 niche-enriched gene clusters revealed that mapped gene clusters (31.2%) were predominantly involved in putative metabolic functions and pathways (Fig. ​(Fig.1e1e and Supplementary Table 3).

Open in a separate window

Extended Data Fig. 1

Fn genetic characterization by niche and subspecies.

a, Size of the Fn pangenome split by the core genome (≥95%) (black) and accessory (<95%) genome (grey). n=10,000 random subsamplings of 135 Fn genomes. Data is plotted as median±s.d. b, Size of the Fn pangenome split by CRC-associated (orange) and oral-associated (blue) niche origin, with respective core and accessory genomes labeled. n=10,000 random subsamplings of 55 Fn CRC-associated and 80 Fn oral-associated genomes. Data is plotted as median±s.d. c, Size of the Fn pangenome split by Fn subspecies, Fna (red), Fnn (gold), Fnp (purple), Fnv (brown), with respective core and accessory genomes labeled. n=10,000 random subsamplings of 75 Fna, 17 Fnn, 33 Fnp, and 10 Fnv genomes. Data is plotted as median±s.d. d-e, Maximum-likelihood dendrograms of d, fadA and e, fplA nucleotide and amino acid sequences. For each genome, tree end points indicate Fn subspecies; Fna (red), Fnn (gold), Fnp (purple), Fnv (brown). f, Column graph depicts the proportion of Fn genomes containing canonical Fn virulence factors, subset by Fna (red), non-Fna (black) subspecies. Statistical analysis performed using two sample Z test, two-tailed. NS, not significant.

Previously published studies using orthogonal approaches have observed a differential distribution of Fn subspecies in tumour tissue from patients with CRC and mucosal biopsy specimens from patients with inflammatory bowel disease^28,29 (IBD). As there is sufficient genetic heterogeneity between the four Fn subspecies that reclassification into separate species has been proposed³⁰, we increased the resolution of our analyses to the subspecies level (Fna, Fn subspecies nucleatum (Fnn), Fn subspecies polymorphum (Fnp) and Fn subspecies vincentii (Fnv); Methods). It is notable that studies seeking to delineate the contributions of Fn in CRC predominately use the model strains Fnn ATCC25586 (ref.⁷), an oral isolate, Fnn ATCC23726 (refs.^20,31), a urogenital isolate, and Fna 7_1 (ref.²⁹), an isolate from a patient with IBD. Here, analysis of the proportion of Fn subspecies by niche found that of the four Fn subspecies, only Fna is significantly associated with the CRC niche (two-sample z-test, two-tailed, P=0.0232; Fig. ​Fig.1f),1f), validating previous studies²⁸, whereas Fnn is significantly enriched in the oral niche (two-sample z-test, two-tailed, P=0.00932; Fig. ​Fig.1f).1f). Thus, we reasoned that the repertoire of genetic factors associated with colonization or virulence in the CRC niche would not be fully represented in Fnn model strains and indeed show that Fnn ATCC25586 contains only 17.60% of CRC-enriched gene clusters (Supplementary Table 3). This supports the use of Fna strains such as Fna 7_1, which has been shown to induce colonic tumours in mouse models³, for mechanistic studies. Comparison of Fn subspecies pangenomes showed that Fna has the smallest core genome compared to those of other subspecies, suggestive of further unresolved Fna genetic heterogeneity (Extended Data Fig. ​Fig.1c1c).

As Fna is enriched in the CRC tumour niche, we tested whether Fn virulence factors previously described as important for host colonization are more prevalent in Fna than in other Fn subspecies. Fn type Va autotransporter virulence factors include fusolisin³², a serine protease that damages host tissue and inactivates immune effectors, FplA (ref.³³), a phospholipase autotransporter that binds to host phosphoinosite-signalling lipids, and the Fap2 (ref.³⁴), Aim1 (ref.³⁵), RadD (ref.³⁶) and CmpA (ref.³⁷) adhesins that mediate Fn interactions with either host cells or other bacterial species. As the role of fplA, aim1, radD and cmpA in CRC remains unclear, to ensure a comprehensive analysis, we queried the presence of these virulence factors across our Fn genome collection. An additional adhesin, FadA (refs.^38,39), mediates Fn attachment to and invasion of host epithelial cells^40,41, with two additional FadA homologues recently identified⁴². Previously analysed fadA distribution in a limited number of Fusobacterium genomes has suggested fadA absence from passively invading species and increased incidence in highly invasive species¹³, although this distribution does not always coincide with in vivo invasion assays^13,38,42,43. We found that although fplA and fadA are well conserved in Fn, their nucleotide and amino acid sequences segregate by Fn subspecies, perhaps indicative of variable interactions with host ligands (Extended Data Fig. 1d,e). Our results show that none of these canonical virulence factors is significantly associated with Fna compared to other Fn subspecies, suggesting that additional unknown genetic factors facilitate the enrichment of Fna in CRC (Extended Data Fig. ​Fig.1f).1f). However, we observed that a subset of Fna strains lacked fap2, cmpA and fusolisin and, by rpoB gene analysis, these Fna strains seemed to form a distinct Fna clade (Fig. ​(Fig.1g1g).

An Fna clade enriched in the CRC niche

To further examine the observation that Fna strains form two distinct clades, we compared phylogenetic trees of housekeeping genes previously used for Fn subspecies typing^28,44. Analysis of these single-marker genes supported the observation that two distinct Fna clades exist, which we call Fna clade 1 (Fna C1) and Fna clade 2 (Fna C2; Extended Data Fig. ​Fig.2a).2a). Beyond these single-marker genes, genome-wide differences between Fna clades are supported by a kSNP⁴⁵ reference-free whole-genome phylogeny (Fig. ​(Fig.2a).2a). To quantify the relatedness of Fna C1 to Fna C2, we compared the average nucleotide identity (ANI), a well-established index that measures the percentage of similarity between genomes, with an established 95% species threshold⁴⁶. Between Fna clades the ANI ranged from 91.61% to 93.11%, comparable to the ANI between other Fn subspecies (Fig. ​(Fig.2b2b and Supplementary Tables 4 and 5). Further, we visualized the patterns of protein-coding genes present across Fna genomes using the Genes in Genomes-Map (GiG-map) tool and found that Fna C1 and Fna C2 have distinct protein-coding gene content (Fig. ​(Fig.2c).2c). This was further supported by principal component analysis (PCA) of Anvi’o gene cluster presence–absence (Fig. ​(Fig.2d).2d). Notably, the frequently used Fna 7_1 strain groups with Fna C2 (Extended Data Fig. ​Fig.2b2b).

Open in a separate window

Fig. 2

Genetic and epigenetic characteristics of Fna clades.

a, A kSNP⁴⁵ maximum-likelihood whole-genome phylogenetic tree. For each Fn genome (n=135), the tree end points indicate the niche origin (CRC (orange); oral (blue)) and the bar colour indicates the Fn subspecies (Fna (red); Fnn (gold); Fnp (purple); Fnv (brown)). Within Fna, the background colour indicates the Fna clade (Fna C1 (green); Fna C2 (lavender)). b, A clustered dendrogram of the ANI matrix. The bar colour indicates the Fn subspecies (Fna (red); Fnn (gold); Fnp (purple); Fnv (brown)). The Fna clades are highlighted with green and lavender boxes. ANI values are reported in Supplementary Tables 4 and 5. c, A GiG-map visualization of the protein-coding gene content across Fn genomes. The top bar colour indicates the Fn subspecies (Fna (red); Fnn (gold); Fnp (purple); Fnv (brown)) and the box colour indicates the Fna clade (Fna C1 (green); Fna C2 (lavender)). The inset on the right highlights groups of protein-coding genes that are distinct between Fna C1 and Fna C2. An interactive GiG-map dataset is available at https://fredhutch.github.io/fusopangea/. d, PCA of Anvi’o gene clusters by presence and absence in each genome. The colours indicate the Fn subspecies and Fna clades (Fnn (gold); Fnp (purple); Fnv (brown); Fna C1 (green); Fna C2 (lavender)). The ellipses are drawn to 95% confidence. e, Left: PCA of Fna genome-wide methyl-modified nucleotide sequences. The ellipses are drawn to 95% confidence. The overlay of the PCA biplot shows the top five nucleotide motifs that are driving the Fna clade bifurcation. Right: a table indicating the distribution of each motif across the Fna clades. The colour indicates the Fna clade (Fna C1 (green); Fna C2 (lavender)). f, A column graph depicting the proportion of Fna CRC-associated and Fna oral-associated genomes, subset by Fna clade (Fna C1 (green); Fna C2 (lavender)). The statistical analysis was carried out using a two-sample z-test, two-tailed. NS, not significant. The graphics in f were created using BioRender.com.

Open in a separate window

Extended Data Fig. 2

Phylogenetics and genomics of Fna clades.

a, Maximum-likelihood phylogenetic trees of Fn single marker genes, 16S rRNA, rpoB, zinc protease, nusA and nusG. For each genome (n=135), tree end points indicate niche origin, CRC (orange) or oral (blue), bar color indicates Fn subspecies (Fna (red), Fnn (gold), Fnp (purple), Fnv (brown)), and background color indicates Fna clades (Fna C1 (green) and Fna C2 (lavender)). b, Genes-in-genomes map (GiG-map) visualization of protein coding gene content across Fna genomes. Boxes highlight Fna clades, Fna C1 (green) and Fna C2 (lavender). Previously published NCBI genomes are labeled by strain name. c, Size of the Fna pangenome split by Fna clade, Fna C1 (green) and Fna C2 (lavender), with respective core and accessory genomes labeled. n=10,000 random subsamplings of 24 Fna C1 and 51 Fna C2 genomes. Data is plotted as median±s.d. d, Composition of Fna pangenome subset by clade. Anvi’o gene cluster (GC) prevalence was used to define core (≥95%), accessory (≥5% and <95%), and rare (<5%) features conserved in both Fna C1 and Fna C2 strains (“Fna core” (≥95% in all Fna strains), “Fna cloud” (≥5% and <95% in all Fna strains), “Fna rare” (<5% in all strains) or unique for strains from each clade. Disparate features are those that do not fall into any of the other noted bins. Plot box shows 25^th percentile, median, and 75^th percentile. Plot whiskers indicate minima and maxima. e, Column graph indicates chromosome sizes in Fna C1 (n=24) and Fna C2 (n=51). Data is plotted as mean±s.e.m. Statistical analysis performed using Welch’s T-test, two-tailed. f, Column graph depicts the proportion of Fna genomes containing innate bacterial genetic defense systems, subset by Fna clades, Fna C1 (green) and Fna C2 (lavender). Statistical analysis performed using two sample Z test, two-tailed. NS, not significant. g, Graph shows the percent relative abundance of Fna C1 (green) and Fna C2 (lavender) in paired saliva (circle) or tumor biopsy (triangle) samples from 39 patients with colorectal adenocarcinomas⁴⁷. Data is plotted as mean±s.e.m. Statistical analysis performed using Welch’s T-test, paired. NS, not significant.

Therefore, we reassessed the genetic, epigenetic and ecological properties of Fna as two genetically distinct clades. Comparison of the Fna clade pangenomes showed that Fna C1 and Fna C2 had similar core genome sizes, although Fna C2 had a larger accessory genome (Extended Data Fig. 2c,d), suggesting that Fna C2 strains harbour additional genetic factors that may be beneficial during colonization of CRC tumours. Consistent with this, comparisons of individual genome size and content indicated that Fna C2 strains have significantly larger chromosome sizes (Welch’s t-test, two-tailed, P<0.00001; Extended Data Fig. ​Fig.2e),2e), more extrachromosomal plasmids (Supplementary Tables 1 and 6) and a greater number of innate genetic defences and mobile genetic elements (Extended Data Fig. ​Fig.2f)2f) compared to Fna C1 strains. PCA analysis of Fna methylomes indicated that Fna C1 and Fna C2 are also epigenetically distinct (Fig. ​(Fig.2e).2e). The methyl-modified DNA motifs most influencing this epigenetic bifurcation are GTN^m6AC (100% Fna C1, 0% Fna C2), GC^m6AG (100% Fna C1, 0% Fna C2) and G^m6ANTC (0% Fna C1, 63% Fna C2; Fig. ​Fig.2e2e and Supplementary Table 7). Although both Fna clades are present in the oral cavity with non-significant differences, only Fna C2 is significantly associated with the CRC niche (two-sample z-test, two-tailed, P<0.00001; Fig. ​Fig.2f).2f). We further tested this observation on publicly available 16S rRNA gene sequencing data from paired tumour tissue and saliva samples from patients with CRC⁴⁷. To gain resolution to the Fna clade level, we identified Fna clade-specific amplicon sequence variants⁴⁸ (Methods and Supplementary Table 8). Supportive of our observations, our data show that Fna C2 is significantly enriched in tumour samples compared to Fna C1 (t-test, paired, P=0.047; Extended Data Fig. ​Fig.2g).2g). However, there was no statistically significant difference between Fna clades in paired oral samples, indicating that although both Fna clades are present in the oral cavity of patients with CRC, only Fna C2 is enriched in the tumour niche (Extended Data Fig. ​Fig.2g2g).

Human lower GI tract Fna C2 enrichment

Pangenome analysis revealed that Fna is composed of two distinct clades but only Fna C2 is enriched in the CRC niche. Low levels of Fna C1 in this niche could be due to poor tumour colonization potential, an inherent lack of virulence or tumour-supportive factors that are possessed by Fna C2, or Fna C1 being unable to evade immune clearance. To interrogate these possibilities and reveal Fna clade-specific genetic factors, we applied a comprehensive inter-Fna clade comparative analysis across all 75 Fna genomes (24 Fna C1, 51 Fna C2). Canonical Fn virulence factors including the adhesins encoded by radD, aim1 and fadA2 were significantly enriched in Fna C1 compared to Fna C2 (two-sample z-test, two-tailed, P<0.00001; Figs. ​Figs.1g1g and 3a) suggesting that their role may be particularly important in the oral cavity. Conversely, as noted, fap2, cmpA and fusolisin are absent from Fna C1 and significantly associated with Fna C2 (two-sample z-test, two-tailed, P<0.00001; Figs. ​Figs.1g1g and 3a). Given the reported epithelial³¹ and immune cell⁴⁹ interactions of Fap2 in CRC, its association with Fna C2 supports its increased adherence and invasion potential in this niche. Co-culture of Fna strains from each clade with a human colon cancer cell line (HCT116) demonstrate that Fna C2 has a significantly higher level of cancer epithelial cell invasion compared to Fna C1 strains (Fig. ​(Fig.3b3b and Extended Data Fig. 3a,b; Welch’s t-test, two-tailed, P=0.0113), indicative of differential invasion potential and/or aerotolerance of individual strains (Extended Data Fig. ​Fig.3c).3c). Further, the Fna clades are morphologically distinct, with Fna C2 cells being significantly longer (Fna C1: 2.01μm average, Fna C2: 5.26μm average) and thinner (Fna C1: 0.39μm average, Fna C2: 0.33μm average) than Fna C1 cells (Extended Data Fig. ​Fig.3d;3d; Welch’s t-test, two-tailed, P<0.00001 length, P<0.00001 width). As bacterial morphology can affect colonization of host niches and susceptibility to host defences⁵⁰, physical differences between Fna clade cells are noteworthy.

Open in a separate window

Fig. 3

Inter Fna-clade comparative analyses.

a, A gene presence versus absence column graph depicting the proportion of Fna genomes containing canonical Fn virulence factors, subset by clade (Fna C1 (green); Fna C2 (lavender)). The statistical analysis was carried out using a two-sample z-test, two-tailed. b, Left: computational confocal analysis of colon cancer epithelial cells (HCT116; grey) co-incubated with representative Fna C1 (green) or Fna C2 (lavender) strains. Scale bars, 4μm. Right: a bar plot demonstrating the percentage of HCT116 cells with intracellular Fna; n=3 biological replicates with 3 analysed z-stacks. Data are plotted as mean±s.e.m. The statistical analysis was carried out using Welch’s t-test, two-tailed. c, A PPanGGOLiN⁵¹ map of the Fna pangenome. Each node represents a gene group, syntenic nodes represent neighbouring genes, the size indicates relative presence across Fna genomes, and the colour depicts pangenome partition (Fna core (red); Fna C1 accessory genome (green); Fna C2 accessory genome (lavender)). The white arrows (left panel) and oblongs (right panel) indicate Fna C2-associated putative eut and pdu operons. d, Schematics of these Fna C2 operons. An interactive PPanGGOLiN map is available at https://fredhutch.github.io/fusopangea/. e,f, Differentially expressed genes (log₂[fold change]≥0.58 and ≤−0.58 with −log₁₀[P value]≥1.30) in a representative Fna C2 strain, SB010, exposed to EA (e) or 1,2-PD (f) compared to unexposed SB010 control. To highlight SB010-unique content, genes also differentially expressed under the same exposure conditions in a representative Fna C1 strain, KCOM 3764, have been removed (Extended Data Fig. 4b,c). The vertical dotted lines indicate the threshold of significant gene expression, defined as log₂[fold change]≥0.58 and ≤−0.58. The statistical analysis was carried out using glmQLFTest, two-sided. The data point colours indicate Fna core (red), Fna C1 accessory (green), Fna C2 accessory (lavender) or Fna cloud (black; present in ≥5% and <95% in all Fna strains) genes. The stars indicate eut and pdu operon genes.

Open in a separate window

Extended Data Fig. 3

Morphological and genomic differences between Fna clades.

Representative Fna C1 and Fna C2 strains co-cultured with human colon cancer cells (HCT116). a-b, Computational analysis of confocal imaging. Independent masks for cancer epithelial cells (grey), and intracellular bacterial cells (Fna C1 green; Fna C2 lavender) were generated. Masks were used to calculate the percent of HCT116 cells with intracellular Fna (Fig. ​(Fig.3b)3b) (see Methods). Scale bar is 20 μm. c, Bacterial aerotolerance was assessed through serial dilution plating at start, mid-point, and endpoint of co-culture. Graph shows resulting bacterial colony forming units per millilter, standardized to start point for each strain. Dashed line indicates normalization equal to one. Statistical analysis performed using a Welch’s T-test, two-tailed. d, Bar plots indicate Fna cell length and cell width as measured from confocal microscopy images, subset by Fna clades, Fna C1 (green) or Fna C2 (lavender); n=45 cells per Fna clade. Data is plotted as mean±s.e.m. Statistical analysis performed using Welch’s T-test, two-tailed. e, KofamKOALA KEGG²⁷ ortholog mapping of Fna clade-enriched gene clusters.

We reasoned that delineation of Fna clade-unique genome content could reveal hitherto unknown genetic factors enabling Fna C2 transit to and survival within the human colonic niche. Predominant genetic differences between Fna clades are consistent with Fna C2 having increased nutrient scavenging mechanisms and enhanced metabolic potential (Extended Data Fig. ​Fig.3e3e and Supplementary Table 9). As functionally related bacterial genes often form co-regulated units (operons), we implemented the Partitioned PanGenome Graph of Linked Neighbors (PPanGGOLiN)⁵¹ tool to assess whether Fna clade-unique genetic factors formed putative operons (Fig. ​(Fig.3c).3c). Consistent with Anvi’o analysis, the PPanGGOLiN analysis showed that Fna C2 syntenic blocks were predominantly associated with metabolic mechanisms (Supplementary Table 10). Thus, the pathoadaptation of Fna C2 to the CRC niche is multifactorial, and in addition to canonical Fn virulence factors is potentially facilitated by enhanced metabolic capabilities.

To validate our inter-Fna clade pangenomic approach, we focused on two Fna C2-associated putative operons consistent with ethanolamine (EA) metabolism (eut) and 1,2-propanediol (1,2-PD) metabolism (pdu; Fig. ​Fig.3d).3d). These operons contribute about 20% of Fna C2-unique gene content (Supplementary Table 10). Enteric pathogens not only gain a competitive advantage through direct metabolism of EA and 1,2-PD, but also exploit their intestinal specificity. Sensing through eut and pdu activates global regulators of virulence and induces transcriptional profiles consistent with GI niche adaptation⁵². Analysis of stool metagenomic datasets from publicly available cohorts of patients with CRC (n=627) and healthy individuals (n=619) indicates that the eut and pdu operons are significantly enriched in patients with CRC (two-sample z-test, two-tailed, eut P<0.00001, pdu P<0.00001; Extended Data Fig. ​Fig.4a4a).

Open in a separate window

Extended Data Fig. 4

Fna clade transcriptomic responses to intestinal metabolites.

a, Bar plots demonstrate the proportion of stool metagenomic samples from patients with CRC or healthy controls in which putative eut and pdu operons were detected. Statistical analysis performed using two sample Z test, two-tailed. Cohort sample sizes are indicated at the bottom of each panel. b-c, Differentially expressed genes (with log2-transformed fold change≥0.58 and ≤ −0.58 with -log10(p-value)≥1.30) of a representative Fna C1 strain (KCOM 3764) and a representative Fna C2 strain (SB010) under (b) ethanolamine (EA) or (c) 1,2-propanediol (1,2-PD) exposure as compared to their respective unexposed controls. Top five significant (-log10(p-value)≥1.30) upregulated and downregulated genes are labeled. d, Differentally expressed genes (with log2-transformed fold change≥0.58 and ≤ −0.58 with -log10(p-value)≥1.30) in SB010 under Vitamin B12 exposure alone as compared to unexposed control. All differentially expressed genes labeled. For b-d, vertical dotted lines indicate the threshold of significant gene expression, defined as log2-transformed fold change≥0.58 and ≤ −0.58. Statistical analysis performed using glmQLFTest, 2-sided. Data point colors indicate whether a gene is categorized as part of the Fna core genome (red), Fna C1-associated accessory genome (green), or Fna C2-associated accessory genome (lavender) by PPanGGOLiN.

Motivated by the conservation of eut and pdu in Fna C2, and their absence in Fna C1, we assessed global transcriptomic responses of Fna cells following exposure to these intestinal-associated metabolites. RNA sequencing of representative Fna C1 and Fna C2 strains after exposure to EA or 1,2-PD indicated that both Fna clades have significant transcriptomic changes (Extended Data Fig. 4b–d and Supplementary Tables 11–15). As Fna C1 is deficient in both eut and pdu, we reasoned that significant transcriptomic changes following exposure to EA or 1,2-PD in Fna C1 would be independent of these operons. Thus, through a subtractive approach, we focused on differentially expressed genes (t-test, two-tailed, P<0.05) in Fna C2 cells that were not differentially expressed (t-test, two-tailed, P>0.05) in Fna C1 cells. Our results demonstrate that in Fna C2, eut and pdu genes are transcriptionally upregulated in response to EA and 1,2-PD, respectively (Fig. 3e,f).

Furthermore, Fna C2 cells exposed to EA or 1,2-PD significantly upregulated 13.02% of Fna C2-associated genes, including canonical Fn virulence factors. Although present in both clades, radD and aim1 are upregulated in the presence of EA in Fna C2 but not Fna C1 cells (Fig. ​(Fig.3e).3e). Virulence factors uniquely present in Fna C2 are additionally upregulated when Fna C2 cells are exposed to EA (cmpA, fusolisin and fap2) or 1,2-PD (fap2; Fig. 3e,f). The upregulation of Fna C2-associated genes and virulence factors known to be important for interactions with human epithelial cells suggests that following their transit to the human GI tract, sensing of these molecules could induce Fna C2 transcriptional profiles consistent with extra-oral niche adaptation.

This also led us to reconsider how Fna C2 might be translocating to extra-oral tumour niches. Previous studies suggest that oral fusobacteria travel to CRC tumours through the bloodstream during transient bacteremia caused by activities such as daily hygiene practices or dental procedures²⁰. Our identification of transcriptionally active eut and pdu operons suggests that direct descent through the GI tract, with subsequent infiltration of CRC tumours through the lumen, may be an additional pathway used by Fna C2. Yet, for GI transit to be a viable route of dissemination, Fna C2 would need to overcome the deleterious effects of extreme acid stress encountered in the stomach (pH1.5–3.5; Extended Data Fig. ​Fig.5a).5a). Assessing the preferential growth pH, we observed that both Fna clades are sensitive to pH below 4.5 (Extended Data Fig. ​Fig.5b).5b). From pH5.5 to 8.5, Fna C2 strains had a significantly higher level of growth activity compared to Fna C1, with maximum growth activity at pH7 (Extended Data Fig. ​Fig.5b).5b). Under basic conditions (pH9.5–10), Fna C1 strains had significantly higher growth activity compared to Fna C2, with maximum growth activity at pH10 (Extended Data Fig. ​Fig.5b).5b). Pangenome analysis also revealed a putative glutamate-dependent acid resistance (GDAR) system conserved across all Fna C2, but absent in Fna C1 (Extended Data Fig. ​Fig.5c5c and Supplementary Table 10). The GDAR system, found in pathogenic and commensal gut bacteria, is one of the most potent acid resistance mechanisms⁵³, with glutamate being the only component necessary for the system to operate at pH≤3 (Extended Data Fig. ​Fig.5d).5d). Using a colorimetric pH change assay, we tested the conversion of glutamine through glutamate into γ-aminobutyric acid in the presence of Fna C1 and Fna C2 and found that the level of this conversion is significantly higher in the presence of Fna C2 strains (Extended Data Fig. 5e,f).

Open in a separate window

Extended Data Fig. 5

Differences in pH preference and acid resistance across Fna clades.

a, Schematic of potential gastrointestinal route from the oral primary niche to the CRC tumor secondary niche. b, Plot indicates growth activity as measured in Biolog PM10 plates for representative Fna C1 (green) and Fna C2 (lavender) strains. Data is plotted as the normalized average across duplicates. Statistical analysis at each pH performed using Welch’s T-test, two-tailed. NG = no growth. c, PPanGGOLiN⁵¹ map of Fna pangenome. Each node represents a gene group, syntenic nodes represent neighboring genes, size indicates relative presence across Fna genomes, and color depicts elements in the Fna C1-associated accessory genome (green), and the Fna C2-associated accessory genome (lavender). White arrow indicates putative glutamate-dependent acid resistance (gdar) operon. d, Schematic depicts mechanism of GDAR acid resistance system. e, Qualitative and f, quantitative measurements of colorimetric assay measuring pH change indicative of conversion of glutamine to glutamate (yellow to green) and conversion of glutamate to γ-aminobutyric acid (GABA) (green to blue) in the presence of representative Fna C1 and Fna C2 strains. n=3 technical triplicates of 3 biological replicates of 3 strains per Fna clade. Data is plotted as mean±s.e.m and statistical analysis performed using an ANOVA. g, Schematic of experiment testing the effects of pH stress by exposure to simulated gastric fluid (SGF) at pH 3 or SGF supplemented with 10mM glutamate at pH 3. Plates show resulting growth for a representative Fna C1 and a representative Fna C2 strain over the course of an hour exposure, as compared to tryptic soy broth (TSB) control at pH ~6.7. h, Bar plots demonstrate the proportion of stool metagenomic samples from patients with CRC (orange) or healthy controls (black) in which a putative gdar operon was detected. Statistical analysis performed using two sample Z test, two-tailed. Cohort sample sizes are indicated at the bottom of each panel. The graphics in a,d,g were created using BioRender.com.

To further mimic effects of pH stress during gastric transit, we exposed Fna clades to simulated gastric fluid at pH3. Both were non-viable after 10min of exposure to simulated gastric fluid. However, in the presence of supplemented glutamate, Fna C2 survived for an extended period (about 60min), which was not observed for Fna C1 lacking GDAR (Extended Data Fig. ​Fig.5g).5g). Analysis of stool metagenomic datasets indicated that gdar operons are significantly enriched in patients with CRC compared to healthy individuals (two-sample z-test, two-tailed, P<0.00001; Extended Data Fig. ​Fig.5h).5h). Thus, in addition to active eut and pdu systems, differences in pH preference and acid resistance mechanisms may contribute to the ability of Fna C2 to access the GI and tumour niches.

Fna C2 affects intestinal tumorigenesis

As Fna C2-enriched gene clusters were predominantly associated with enhanced metabolic potential (Extended Data Fig. ​Fig.3e3e and Supplementary Table 10), we sought to determine whether Fna treatment of the dextran sodium sulfate-induced colitis Apc^Min+/− mouse model⁵⁴ of CRC affected intestinal tumorigenesis and metabolic pathways in vivo (Fig. ​(Fig.4a).4a). To capture a higher proportion of Fna clade-specific accessory genes (Extended Data Fig. ​Fig.2d),2d), a mix of three representative strains for each clade was used. Following the administration of a single oral gavage of Fna C1 mix, Fna C2 mix or vehicle control, we observed a significant increase in the number of intestinal adenomas in Fna C2-treated mice compared to both Fna C1 and vehicle control independently (Extended Data Fig. 6a,b; analysis of variance (ANOVA), P=0.0065 and P=0.0069, respectively), specifically in the large intestine (Fig. ​(Fig.4b;4b; ANOVA, P=0.0070 and P=0.0009, respectively; Extended Data Fig. ​Fig.6c).6c). There was no significant difference in adenoma burden between Fna C1 treatment and vehicle control mice. Low-level Fn was inconsistently detected during the course of the study (Extended Data Fig. ​Fig.6d).6d). We carried out liquid chromatography–mass spectrometry global metabolomics on intestinal tissue from each treatment arm for comparative metabolite analysis (Supplementary Table 16). Partial least squares discriminant analysis of measured intestinal metabolites demonstrated that Fna C2-treated mice formed a distinct cluster away from other treatment arms, suggesting a differential metabolic profile. However, intestinal metabolites from Fna C1-treated and vehicle control mice had more similar metabolite profiles, clustering together (Fig. ​(Fig.4c4c).

Open in a separate window

Fig. 4

Fna C2 impact on intestinal tumorigenesis and metabolism.

a, A schematic of the study with 6–8-week-old Apc^Min+/− mice receiving streptomycin and dextran sodium sulfate (DSS) treatment to alter the native microbiome and induce colitis, respectively. Mice were orally gavaged with vehicle control (arm 1) or a mix of three representative Fna C1 (arm 2) and Fna C2 (arm 3) strains. A strain mix was used to capture a higher proportion of Fna clade-specific accessory genes (Extended Data Fig. ​Fig.2d).2d). The mice were monitored until the end-point at 6 weeks post-gavage when they reached 15–17 weeks of age. b, A plot indicating the number of adenomas in the large intestine by treatment arm (vehicle control (grey); Fna C1 treated (green); Fna C2 treated (lavender); n=8 mice per arm). The data are plotted as mean±s.e.m. The statistical analysis was carried out using one-way ANOVA. c, Partial least squares discriminant analysis of detected intestinal metabolites (n=1,296). The colours represent the treatment arm (vehicle control (grey); Fna C1 treated (green); Fna C2 treated (lavender)). The graphics in a were created using BioRender.com.

Source Data

Open in a separate window

Extended Data Fig. 6

Intestinal adenoma burden and fecal Fusobacterium load in Fna treated mice.

a, Schematic of study with Apc^Min+/− mice orally gavaged with vehicle control (Arm 1) or representative Fna C1 (Arm 2) and Fna C2 (Arm 3) strains post streptomycin and dextran sodium sulfate (DSS) to clear the native microbiome and induce colitis, respectively. b-c, Plots indicate the number of adenomas in b, the small and large intestines combined and c, the small intestine by treatment arm, vehicle control (grey), Fna C1-treated (green), or Fna C2-treated (lavender). Data is plotted as mean±s.e.m; n=8 mice per treatment arm. Statistical analysis performed using one-way ANOVA. d, Fusobacterium-targeted qPCR was performed on fecal pellets from one-day post-gavage (PG) to study endpoint. Fn presence in each sample is plotted as Fn copies per ng of fecal DNA with a detection limit (D.L.) of 1 Fn copy. Each sample was run in triplicate, and samples with ≥2 duplicates with detectable signal are included. The graphics in a were created using BioRender.com.

Source Data

Of 1,296 metabolites measured (Extended Data Fig. 7a–c and Supplementary Table 16), comparative analysis demonstrated a significant enrichment in glutathione metabolism and γ-glutamyl amino acid pathways in Fna C2-treated mice, compared to both Fna C1-treated and vehicle control mice (Extended Data Figs. 7a–c, ​,8a8a and ​and99 and Supplementary Table 17). Specifically, we observed a significant increase in the levels of precursors to γ-glutamyl-cysteinyl-glycine (GSH) synthesis, including cysteine and γ-glutamylcysteine, decreased levels of glutathione in its reduced form (GSH) and significantly higher levels of the GSH degradation product 5-oxoproline (Extended Data Figs. ​Figs.8b8b and ​and99 and Supplementary Table 16). Consistent with γ-glutamyl amino acid generation from reduced glutathione in the presence of γ-glutamyl transpeptidase, increased levels of γ-glutamyl amino acids and cysteine, glycine and cysteinyl-glycine were also observed (Extended Data Fig. 8b,c and Supplementary Table 16). GSH deficiency or an elevated ratio of oxidized (GSSG) to reduced (GSH) forms of glutathione increases the vulnerability of mammalian cells to oxidative stress, inflammation and tumour progression⁵⁵. The GSSG/GSH ratio of Fna C2-treated mice significantly increased by 3.5- and 3.0-fold compared to the control and Fna C1-treated groups respectively, suggesting increased oxidative stress (Extended Data Fig. ​Fig.8d;8d; one-way ANOVA, P=0.0031 and P=0.0047, respectively). Studies have demonstrated that metabolism of GSH by γ-glutamyl transpeptidase can exert pro-oxidant effects. In cancer cells, this is a source of endogenous reactive oxygen species that can facilitate persistent oxidative stress and contribute to genomic instability^56,57. Consistent with this, significantly increased levels of other markers of oxidative stress including cystine and cysteine-glutathione disulfide and significantly decreased levels of polyamines capable of scavenging reactive oxygen species, including putrescine, spermidine and spermine, were observed in Fna C2-treated mice (Extended Data Figs. ​Figs.8b8b and 9). Notably, recent work from our group demonstrated that Fusobacterium was predominantly associated with epithelial cells harbouring severe chromosomal abnormalities¹⁰, one of the most common forms of genomic instability in cancer.

Open in a separate window

Extended Data Fig. 7

Intestinal metabolite changes in Fna treated mice.

Scatter plots shows fold change for 1,296 detected metabolites from a, Fna C2 treated mice compared to control, b, Fna C2 treated mice compared to Fna C1 treated mice, and c, Fna C1 treated mice compared to control colored by metabolic categorization. Top five characterized metabolites are labeled. Plots below indicate the pathway enrichment score all enriched pathways. Bars are colored by metabolic categorization.

Open in a separate window

Extended Data Fig. 8

Altered metabolites in Fna treated mice.

a, Plots show the log2-transformed fold change in metabolite level between treatment arms by the -log10(p-value). Points indicate individual metabolites colored by their metabolic categorization. Dotted lines indicate the threshold of significant gene expression, defined as log2-transformed fold change≥0.58 and ≤ −0.58 (vertical lines) and a -log10(p-value)≥1.30 (horizontal line). Statistical analysis performed on natural log-transformed values using one-sided T-test. b, Clustered heatmap of top fifty metabolites across study arms. Dendrogram groups individual samples by similarity of metabolite profile. Bar color for each sample indicates treatment arm, vehicle control (grey), Fna C1-treated (green), or Fna C2-treated (lavender). c, Schematic of glutathione metabolic pathway. d, Plot demonstrates the ratio between oxidized (GSSG) and reduced (GSH) glutathione levels for each treatment arm, vehicle control (grey), Fna C1-treated (green), or Fna C2-treated (lavender); n=4 mice per arm. Data is plotted as mean±s.e.m. Statistical analysis performed using one-way ANOVA. e, Schematic of eicosanoid metabolic pathway.

Open in a separate window

Extended Data Fig. 9

Significantly altered intestinal metabolites in Fna C2-treated mice.

Plots of individual, characterized metabolites, with metabolic categories having > 10 significantly altered metabolites shown. Each plot shows the log2-transformed log change between Fna C2-treated mice compared to Fna C1-treated mice versus the log2-transformed log change between Fna C2-treated mice compared to control mice. Dashed lines indicate the threshold of significantly altered metabolites, defined as log2-transformed fold change (FC)≥0.58 and ≤ −0.58. Upper right quadrant indicates metabolites that are significantly elevated in Fna C2-treated mice as compared to both Fna C1-treated and control mice. Lower left quadrant indicates metabolites that are significantly lower in Fna C2-treated mice as compared to both Fna C1-treated and control mice. Top metabolites in each of these quadrants are labeled.

In addition to their role in combating oxidative stress, polyamines are able to suppress inflammation through inhibition of macrophage cytokine synthesis⁵⁸. Consistent with increased inflammation, significantly higher levels of N-monomethylarginine and dimethylarginine were observed in Fna C2-treated mice compared to other treatment arms (Extended Data Fig. ​Fig.8b).8b). Both of these metabolites inhibit the synthesis of the anti-inflammatory agent nitric oxide. Furthermore, we observed significantly higher levels of pro-inflammatory prostaglandins and ceramides, including prostaglandin A2, N-palmitoyl-sphingosine and N-palmitoyl-sphingadienine (Extended Data Fig. ​Fig.8b).8b). Ceramides can also be metabolized by cancer cells to reduce tumour cell apoptosis and proliferation⁵⁹. Other metabolites that promote cancer cell proliferation and metastasis across a range of cancers include eicosanoids, which are similarly significantly increased in Fna C2-treated mice compared to Fna C1-treated or vehicle control mice (Extended Data Figs. 7a–c and 8e and Supplementary Table 17). This included increased levels of 6-keto prostaglandin F1-α (Extended Data Fig. ​Fig.8b)8b) through COX2 (also known as PTGS2) metabolism of arachidonic acid. Notably, COX2 was previously reported to be one of the most upregulated genes in Fn-associated human colorectal tumours³. Overall, our results demonstrate the ability of Fna C2, but not Fna C1, to metabolically affect the intestinal milieu towards pro-oncogenic conditions.

Fna C2 enrichment in human CRC cohorts

As Fna C2 strains are both significantly enriched in the CRC niche (Fig. ​(Fig.2f)2f) and increase intestinal tumorigenesis in our mouse model compared to Fna C1 (Fig. ​(Fig.4b4b and Extended Data Fig. 6a–c), we next sought to determine the prevalence and abundance of these Fna clades in human tissue and stool specimens through culture-independent approaches. We carried out bacterial 16S rRNA gene sequencing on resected tumour tissue from 116 patients with treatment-naive CRC (CRC cohort 1) and on adjacent normal tissue from 62 of these patients (Supplementary Table 18). Comparing the percentage relative abundance of different Fusobacterium species between paired tumour and adjacent normal tissue (n=62 patients), we observed that Fn was the only Fusobacterium species significantly enriched in tumour tissue compared to adjacent normal (Fig. ​(Fig.5a5a and Supplementary Table 19), supportive of previous reports^2,60 (t-test, paired, P=0.0022). However, using the Fna clade-specific amplicon sequence variants to resolve Fn to a higher taxonomic resolution that includes Fna C1, Fna C2 and non-Fna subspecies of Fn, we demonstrate that only Fna C2 is significantly enriched in tumour compared to paired normal tissue (Fig. ​(Fig.5a,5a, and Supplementary Tables 8 and 19; t-test, paired, P=0.0093). As neither Fna C1 nor non-Fna subspecies of Fn are significantly enriched, this suggests that it is specifically Fna C2 that is driving the previously reported enrichment of Fn in human CRC tumours. Furthermore, across two independent patient cohorts (CRC cohort 1 n=116 and CRC cohort 2 n=86), we demonstrate that within CRC tumour tissue, Fna C2 is significantly enriched compared to Fna C1 (Fig. ​(Fig.5b5b and Supplementary Tables 20 and 21; t-test, paired, cohort 1 P=0.0009, cohort 2 P=0.0014), supporting our observations at the Fna strain level (Fig. ​(Fig.2f2f).

Open in a separate window

Fig. 5

Fn in human tissue microbiome and stool metagenomic specimens.

a, Plots showing the relative abundance for Fusobacterium species (Fg, F.gonidiaformans; Fh, F.hwasookii; Fm, F.mortiferum; Fnavi, F.naviforme; left plot), and Fn subspecies and Fna clades (right plot) using microbial 16S rRNA gene sequencing of paired tumour (orange) and normal adjacent (black) tissue (n=62 patients with CRC). Amplicon sequence variants were used to obtain Fna clade resolution (Extended Data Fig. ​Fig.1010 and Supplementary Table 8). The data are plotted as mean±s.e.m. The statistical analysis was carried out using one-sided t-test, paired. b, Plots showing the relative abundance for Fna C1 (green) and Fna C2 (lavender) within patient primary colorectal tumour tissue from two independent cohorts (cohort 1 (n=116) this study; cohort 2 (n=86) BioProject PRJNA362951). The data are plotted as mean±s.e.m. The statistical analysis was carried out using one-sided t-test, paired. c, Fna C1 and Fna C2 detection in stool metagenomic data from patients with CRC and healthy individuals. The left plot shows the pooled effect sizes for Fna C1 (green) and Fna C2 (lavender) calculated using a meta-analysis of standardized mean differences and a random-effects model on MetaPhlAn4 (ref.⁶³) species-level genome bin abundances on all CRC samples (n=627) and samples from healthy individual (n=619). The right plot shows the effect sizes for Fna C1 and Fna C2 calculated using the same approach, but here samples in which Fna C1 co-occurred with Fna C2 were excluded. The data are plotted as mean±s.e.m. The statistical significance was assessed by a Wald test, two-sided. All P values are corrected using the Benjamini–Yakuteli method. d, Fna C1 and Fna C2 presence in stool metagenomes of patients with CRC. The bars indicate individual stool samples from patients with CRC (n=627) and are coloured by Fna C1 and Fna C2 detection (Fna C1 detected (green); Fna C2 detected (lavender); Fna undetected (grey)). The lower brackets indicate the number of stool samples in which Fna C1 occurred independently (n=5), Fna C2 occurred independently (n=147), Fna clades co-occurred (n=31) or Fna clades were not detected (n=444). The graphics in a–c were created using BioRender.com.

We next sought to determine whether the prevalence of Fna clades differed between patients with CRC and healthy individuals. To do so, we analysed stool metagenomic datasets from publicly available cohorts of patients with CRC (n=627) and healthy individuals (n=619; Extended Data Fig. ​Fig.1010 and Supplementary Table 22). Fna was detected in 29.2% of stool samples from patients with CRC and 4.8% of stool samples from healthy individuals (Supplementary Table 23). Meta-analysis of standardized mean differences by a random-effects model for Fna C1 and Fna C2 demonstrated that both Fna clades have a significant pooled effect size associated with CRC (Fna C1 effect size=0.21, 95% confidence interval (0.09,0.32), P=4.45×10⁻⁴; Fna C2 effect size=0.45, 95% confidence interval (0.34,0.56), P=5.55×10⁻¹⁵; Fig. ​Fig.5c,5c, Extended Data Fig. ​Fig.1111 and Supplementary Table 24). However, the effect size for Fna C2 was larger than that for Fna C1. Notably, in the absence of Fna C2 co-occurrence, Fna C1 was not significantly associated with CRC (Fig. ​(Fig.5c,5c, Extended Data Fig. ​Fig.1111 and Supplementary Table 25). Although synergistic interactions between CRC-enriched microbes have previously been reported⁶¹, it is not clear whether Fna C1 co-occurrence with Fna C2 results in a compounding pathogenic effect. Similar to our observation in CRC tumour tissue (Supplementary Tables 20 and 21), our data show that Fna C2 is more prevalent and abundant in the stool of patients with CRC than Fna C1 (Fig. ​(Fig.5d5d and Extended Data Fig. ​Fig.10)10) and is furthermore the only Fn subgroup significantly enriched in the stool of patients with CRC compared to healthy individuals (Extended Data Fig. 10a). These culture-independent human specimen analyses support our strain-level genomic discovery that Fna C2 is the dominant CRC-associated Fna clade (Fig. ​(Fig.2).2). This further highlights the significance of our in vitro (Fig. ​(Fig.3)3) and in vivo (Fig. ​(Fig.4)4) findings demonstrating the increased virulence and tumorigenic potential of Fna C2 compared to Fna C1.

Open in a separate window

Extended Data Fig. 10

Detection of Fn subspecies and Fna clades in human stool metagenomes.

Detection of a, Fn subspecies (Fnn, Fnv, Fnp) and a-k, Fna clades (Fna C1 and Fna C2) in stool metagenomic data from previously published independent cohorts from patients with CRC and healthy controls. a, Plot shows the percent relative abundance of each Fn subspecies and Fna clade in stool samples from patients with CRC (n=627) or healthy controls (n=619). Data is plotted as mean±s.e.m and statistical analysis performed using a one-way ANOVA. b-k, Samples are plotted both by (b-k) individual cohorts and (k) pooled results. In each panel, top plot shows the percent relative abundances of Fna C1 (green) and Fna C2 (lavender) in each sample. Data is plotted as mean±s.e.m and statistical analysis performed using Welch’s T-test, paired. Bottom plot demonstrates the proportion of stool samples from patients with CRC and the proportion of stool samples from healthy controls in which Fna C1 and Fna C2 were detected. Data is plotted as mean + s.d. and statistical analysis performed using two sample Z test, two-tailed. Cohort sample sizes are indicated at the bottom of each panel.

Open in a separate window

Extended Data Fig. 11

Meta-analysis of Fna in human stool metagenomes.

Fna C1 and Fna C2 detection in stool metagenomic data from previously published independent cohorts from patients with CRC and healthy controls plotted by individual cohorts. Left plot shows the effect sizes for Fna C1 and Fna C2 calculated across all samples (CRC n=627, healthy control n=619) using a meta-analysis of standardized mean differences and a random effects model on MetaPhlAn4⁶³ species-level genome bins (SGB) abundances. Right plot shows the effect sizes for Fna C1 and Fna C2 calculated using the same approach but with samples where Fna C1 co-occurred with Fna C2 excluded (CRC n=596, healthy control n=616). Data is plotted as mean±s.e.m. Statistical significance assessed by Wald test, two-sided. All p-values are corrected via the Benjamini-Yakuteli method.

Discussion

Advances in next-generation sequencing have revealed the presence of bacterial communities within human tumour tissues. A key challenge for cancer microbiome research is to move beyond the characterization of microbial composition in tumours towards functional studies that determine whether, and how, these microbes are contributing to disease. In CRC, Fn gained early and continued attention owing to the fact that this bacterium was rarely detected in the lower GI tract of healthy individuals¹, yet enriched within the CRC tumour microbiome^2,60. Fn species are normal members of the human oral microbiota, and strains from the oral cavity are thought to seed CRC tumours^19,20. However, the noted genetic and phenotypic heterogeneity^12–15 of Fn led to an open question of whether Fn strains that colonize and dominate human tumours harbour distinct genetic attributes that contribute to CRC initiation or progression. Through large-scale culturing, sequencing and comparative genomic analyses of human CRC and non-cancer oral Fn strains, we revealed the distinct CRC-enriched genetic factors of Fn. Further, we identified that these CRC-enriched factors were predominantly present within a specific clade of Fna. This was mirrored by our discovery that Fna is bifurcated into two distinct clades: Fna C1, which is largely restricted to the oral cavity, and Fna C2, which dominates the human CRC tumour niche. Notably, only Fna C2 induced tumours and altered intestinal metabolism towards increased oxidative stress within a CRC animal model. Further, comparative genomic analysis between Fna clades revealed the genetic elements that cumulatively engender the pathoadaptation of Fna C2 to the CRC niche. Given the power of using Fna C1 as a comparative group for Fna C2, we created an interactive website to enable the exploration of Fna pangenomic datasets designated The Fusobacterium Pangenome Atlas at https://fredhutch.github.io/fusopangea/. Collectively, this work demonstrates that Fna C2 is a highly virulent subgroup of Fn that should be the primary focus for mechanistic studies and therapeutic drug design in CRC.

Methods

Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-024-07182-w.

Acknowledgements

Author contributions

M.Z.-R., S.B. and C.D.J. designed the study and wrote the paper. S.B., E.S., K.D.L. and A.G.K. processed patient tissue specimens. M.Z.-R., K.D.L., S.C.Y., S.-N.P., Y.K.L., J.-K.K. and S.B. carried out microbial isolation and culture. M.Z.-R., D.S.J., S.L.C. and C.D.J. carried out DNA isolation, genome sequencing and genome assembly. M.Z.-R., S.S.M., H.B., H.W. and D.S.J. carried out pangenomic analyses on genomes and methylomes. A.B.M. and P.M. analysed stool metagenomes and carried out meta-analysis. M.Z.-R. carried out RNA sequencing experiments and analysis. K.D.L. carried out co-culture assays and subsequent microscopy-based cell invasion analysis. M.Z.-R., K.D.L., Y.W. and A.G.K. carried out mouse trials and downstream processing of mouse samples. M.Z.-R. analysed metabolite data. E.F.M. carried out the glutaminase assays and phenotyping using Biolog PM10 plates. M.Z.-R., A.D.W. and S.B. carried out statistical analyses. M.Z.-R., F.E.D., N.S., S.B. and C.D.J. obtained funding and supervised computational and wet lab experiments. All authors read and provided edits to the paper and contributed to the final version.

Peer review

Data availability

All genomes from this study are available at the NCBI under the BioProject accession number PRJNA549513 and all methylomes are available in the Restriction Enzyme Database (REBASE). Raw sequencing data from RNA-sequencing experiments are available in the NCBI Sequence Read Archive repository under the BioProject accession number PRJNA937266. Raw sequencing data from 16S rRNA sequencing experiments are available in the NCBI Sequence Read Archive repository under the BioProject accession number PRJNA1064180. Source data are provided with this paper.

Competing interests

S.B. has consulted for GlaxoSmithKline and BiomX. C.D.J. has consulted for Series Therapeutics and Azitra. S.B. is an inventor on US patent application no. PCT/US2018/042966, submitted by the Broad Institute and Dana-Farber Cancer Institute, which covers the targeting of Fusobacterium for the treatment of CRC. S.B., C.D.J. and M.Z.-R. are inventors on US patent application no. F053-0188USP1/22-158-US-PSP, submitted by the Fred Hutchinson Cancer Center, which covers the modulation of cancer-associated microbes. K.D.L. is employed by NanoString Technologies at present. The remaining authors declare no competing interests.

Footnotes

Contributor Information

Susan Bullman, Email: gro.hctuhderf@namllubs.

Christopher D. Johnston, Email: gro.hctuhderf@notsnhoj.

References