Evolution of Invasion in a Diverse Set of Fusobacterium Species

Phylogeny and genome content for a highly diverse set of Fusobacterium strains.We generated high-quality draft genome sequences for 21 Fusobacterium strains, representing 7 species. Strains were isolated from a variety of human habitats that were either healthy or inflamed at the time of sample collection (see Table S1 in the supplemental material) (5, 22, 27–29). Additionally, we generated finished genomes for five F. nucleatum strains (see Table S1) that, together with the two previously finished genomes (30, 31), resulted in the inclusion in our study of seven finished F. nucleatum genomes, representing four subspecies.

At the time of sequencing, only partial 16S sequence information was available for 14 isolates, so definitive taxonomic classification was not possible. Using our newly sequenced genomes and 5 previously sequenced Fusobacterium genomes (30–32; http://www.hgsc.bcm.tmc.edu/; http://genome.wustl.edu) (see Table S1 in the supplemental material), we constructed a phylogenetic tree from alignment of 498 orthologous genes, or orthogroups, conserved across all strains (Fig. 1; see Materials and Methods and Table S2 in the supplemental material) and taxonomically classified previously unnamed clinical isolates (see Fig. S1 in the supplemental material). The resolved taxonomy, showing that species fall into three main lineages (Fig. 1), was in general agreement with the taxonomy based on 16S sequences from the Living Tree Project (2). However, our results allowed for a more highly resolved view of evolution within this genus and revealed that there was an adaptive radiation, where the last common ancestor of all Fusobacterium species diversified into three major lineages (see Text S1 in the supplemental material) at an early point in its evolution from Leptotrichia.

• Open in new tab
• Download powerpoint

FIG 1 

Phylogenetic tree based on nucleotide sequences of 498 core orthogroups, or orthogroups containing exactly one copy from each of the 26 Fusobacterium strains plus the outgroup, Leptotrichia buccalis. Bootstrap values are indicated for each node. The node indicated with an arrow illustrates a 3-way trifurcation (based on bootstrap values and individual orthogroup trees [see Text S1 in the supplemental material]), representing an adaptive radiation. The five clades are outlined with boxes. The clades containing species believed to actively or passively invade host cells are indicated with dark and light gray shading, respectively. The mechanism by which F. mortiferum (in white) can invade cells is unknown.

We observed substantial variation in genome sizes and gene content among the sequenced species (see Table S1 and Fig. S2 in the supplemental material). Genome sizes differed by as much as 2 Mb (1.7 to 3.7 Mb). Only 556 orthogroups (24%) were shared among all strains. Members of the same species and subspecies exhibited remarkable gene content plasticity: F. nucleatum strains shared only 59% of their orthogroups, while members of the closely related F. nucleatum subsp. animalis subspecies shared only 70%, despite having >99% nucleotide identity among shared genes. In addition, few orthogroups uniquely defined each species and subspecies. (For example, 1%, or 36 orthogroups, were found exclusively in all 14 F. nucleatum genomes, and 1%, or 32 orthogroups, were found exclusively in all 6 F. nucleatum subsp. animalis genomes.) Average nucleotide identity (ANI) plots based on whole-genome data (33, 34) (see Fig. S1 in the supplemental material) suggested that each F. nucleatum subspecies could be considered a separate species (89 to 93% ANI), and F. periodonticum could also be subdivided into separate species, with some strain comparisons having only 92 to 94% ANI (see Text S1 in the supplemental material), likely contributing to the low numbers of species-defining genes. Species with fewer genomes had more defining orthogroups (e.g., 6%, or 220, were found exclusively in F. ulcerans), likely due to their being underrepresented in our data set (see Table S3 in the supplemental material) and their higher ANI (>95%) (see Text S1). Species-specific orthogroups included genes encoding extracellular features, such as adhesins, membrane-associated transporters, receptors, and extracellular solute-binding proteins (see Text S1). That some adhesins and membrane-associated proteins were species specific has been previously documented in Fusobacterium (35) and points to surface attachment as a driver of the evolution of this genus.

Since the majority of species- and subspecies-defining orthogroups encoded “hypothetical proteins” of unknown function, we compared the sequences of genes against those in functional databases, including the KEGG (36), Pfam (37), and Gene Ontology (GO) (38) databases, to increase our power to assign function to hypothetical proteins. We then searched for broader functional categories expanded in individual species and subspecies that were not obvious by gene annotation and ortholog clustering (see Table S4 and Text S1 in the supplemental material) but that might help to explain what drove speciation and adaptation of Fusobacterium species to different environments. Again, the results indicated that the majority of species-specific orthogroups encoded extracellular or membrane-related proteins, including 46 of the 72 species-specific orthogroups present in either F. nucleatum or F. periodonticum (64%). Expansions were also observed among genes related to amino acid metabolism and cofactor biosynthesis (F. nucleatum subsp. polymorphum and F. periodonticum), gene regulation and signaling (F. ulcerans), YadA-related autotransporter adhesins (F. necrophorum), and species-specific hemagglutinin-related genes (F. nucleatum).

Identification of new determinants associated with active invasion.Only some Fusobacterium species can invade host cells independently (active invaders), while others require either helper organisms, coinfection with viruses, or compromised mucosal integrity (passive invaders) (4–6). The FadA adhesin was previously identified as an important factor in active invasion (7, 8) and has been identified in several actively invading Fusobacterium species (39). To identify other genetic factors associated with active invasion, we first stratified species into active and passive invader clades using a Bayesian approach (see Materials and Methods) and then correlated this structure with information from the literature regarding each species’ invasion potential (4–6, 21). Five distinct clades were identified among Fusobacterium species (Fig. 1). Three clades—clades A (F. nucleatum) and B (F. periodonticum) belonging to lineage 1 and clade C (F. ulcerans/F. varium) belonging to lineage 2—consisted entirely of species known to actively invade host cells (active invaders). Clade D, which included the remaining species in lineage 2, F. mortiferum, was of unknown invasion potential. Clade E constitutes lineage 3 and contained the known passive invader species F. necrophorum and F. gonidiaformans.

Active invaders had genomes that were, on average, 560 kb larger than those of passive invaders and contained 257 more genes (see Table S1 in the supplemental material). A large fraction of the additional genes encoded membrane-related proteins (as defined by GO criteria), doubling the membrane-related protein coding capacity in active invaders. Active invaders also had 1.6-fold more genes with predicted signal peptides (365 versus 233 per genome), suggesting an extracellular role for expanded gene families in active invasion.

When we compared orthogroups between active and passive invaders, 44 were present in all active invaders and absent in all passive invaders (see Table S2b and Table S5 in the supplemental material). Of the genes exclusive to active invaders, 32 were annotated as “hypothetical proteins.” However, we were able to gain clues to the function of many by examining Pfam or GO functional annotations. This highlighted several genes with known links to virulence and pathogenesis, including those encoding branched-chain amino acid transport (40, 41), components of type IV pili (42) and the related bacterial type II secretion system (43–46), a patatin-like phospholipase (47), and META domain-containing (48) and chorismate mutase domain-containing (49, 50) proteins. In addition, 43% (19 of 44) of the active invader-specific genes were predicted to contain either a signal peptide or transmembrane domain, highlighting the importance of evolved surface features on the invasive strains.

To ensure that the maximum number of genes were accounted for in the comparative analysis of these high-quality draft genomes, we used the GO, Pfam, and KEGG functional annotations to identify additional differences in actively versus passively invading species. Using this approach, additional GO-determined functional families were found that had significantly larger numbers of genes in the active invaders, with most being associated with the membrane (Q < 5e−10) (Fig. 2A; see Table S5 in the supplemental material), consistent with our initial analysis identifying a 2-fold enrichment in membrane proteins among this group. Of the seven significantly enriched Pfam domains (Fig. 2B and C), five were associated with adherence, including the FadA adhesin domain (Q < 4.5e−21). By this analysis, FadA family-encoding genes were found to be exclusive to active invader genomes, including both strains of F. ulcerans. It had been suggested previously (39) that FadA was not present in F. ulcerans ATCC 49185; however, we found seven related FadA family genes in that strain and eight in F. ulcerans 12-1B. These were not clustered within the same orthogroup as canonical FadA (FN0264 in F. nucleatum) because they shared only 30 to 38% amino acid identity. The passive invaders, F. necrophorum and F. gonidiaformans, completely lacked FadA family genes, as did F. mortiferum.

• Open in new tab
• Download powerpoint

FIG 2 

Gene categories expanded in active and passive invaders (Q < 0.0005). (A) The top 12 GO terms expanded in the active invaders were largely membrane related. (B) Pfam domains expanded in the active invaders include the known virulence-related adhesins FadA and RadD, as well as a massive expansion of MORN2 domains of unknown function. (C) A different set of adhesins, including the trimeric autotransporter adhesion protein YadA, is expanded in the passive invaders. (See Table S5 in the supplemental material for a full listing of categories expanded in both active and passive invaders for Q < 0.05.) (D) Clade-specific summary of expanded Pfam domains. The number of dots represents the average number of domains present in a genome from each of the five clades. For MORN2, there can be multiple domains per gene.

Interestingly, active invaders also contained 6 times as many genes encoding the adhesin-related autotransporter β-domain (Q < 4.4e−7 [11.0 versus 1.6 genes]) (51), a component of the RadD family of outer membrane protein adhesins known to be involved in pathogenesis of F. nucleatum through cell death in human lymphocytes (51, 52). Bacterial microcompartment (BMC) domains were also found exclusively in the active invaders (Q < 9.9e−5 [4.9 genes]). BMC-containing proteins are associated with polyhedral microcompartment organelles, such as carboxysomes, which give bacteria the ability to adapt to new niches via metabolic innovation and have been previously implicated in pathogenesis (53) (Fig. 2B). While BMC domains have not been directly implicated in adhesion, these genes in F. nucleatum subsp. nucleatum ATCC 25586 were recently shown to be highly induced under conditions that promote aggregation (11), suggesting that these proteins are also involved in aggregation and adhesion.

Of perhaps greatest importance, MORN2 domains (Pfam identifier PF07661) were found to be the most extensively enriched domain among active invaders (Fig. 2B) (Q < 2.5e−21 [32.4 versus 4.8 genes]). None of the 697 genes encoding MORN2 domains (see Table S2c in the supplemental material) within our data set had known function, although 87% were predicted by SignalP (54) to have signal peptides targeting them for export into the extracellular environment or membrane insertion, highlighting a potential role for these domains at the host-pathogen interface. To further validate these results, we verified that all the domains enriched in the active invaders were also enriched in six more recently sequenced F. nucleatum genomes isolated from cancerous tumors (see Text S1 in the supplemental material). We also verified that all 44 of the orthogroups exclusive to active invader strains were also present (see Text S1).

Enriched features in passive invaders included two types of YadA domain proteins (Q < 1.0e−5 [9.6 versus 3.4 copies] and Q < 2.0e−6 [9.7 versus 3.1 copies]). These proteins are found in trimeric autotransporter adhesins implicated in host cell adherence (55). Passive invaders also exclusively encoded proteins with the domain DUF2147 (domain of unknown function) (Q <1.1e−4 [2 copies per species]). Genes with DUF2147 also had predicted signal peptides and were located near other genes that encoded predicted membrane proteins, suggesting roles at the host-microbe interface.

Evolution of active invasion.With few exceptions (see Table S5 in the supplemental material), orthogroups unique to active invading strains were scattered around the genome, with little evidence of having been acquired together. However, genes containing predicted Pfam structural domains enriched among active invaders tended to cluster together (Fig. 3; see Text S1 in the supplemental material), suggesting that they evolved together for related function, such as novel surface feature assembly. This observation was supported by neighborhood analysis, which showed a statistically significant association among orthogroups containing the expanded Pfam families—MORN2, FadA, and RadD—occurring near each other and also near genes encoding other virulence-related proteins and membrane proteins (see Table S6 and Text S1 in the supplemental material). While IS elements and proteins related to phage and transposition were found near regions containing Pfam families expanded in active invaders (see Table S6 and Text S1), there was little evidence of them being part of recognizable prophages or recently acquired through lateral transfer (see Materials and Methods), although one gene containing MORN2 domains was present on a predicted plasmid in F. nucleatum 3_1_27 (see Table S1 in the supplemental material).

• Open in new tab
• Download powerpoint

FIG 3 

Gene families expanded in active invaders clustered in the genome. This multiple alignment of our set of seven finished F. nucleatum genomes shows an example of a region with close physical association among FadA, RadD, and MORN2 family proteins. The gray-shaded connectors represent orthologous protein pairs. This genomic region also contains an ompA gene shown to be involved in biofilm formation in Fusobacterium (69), as well as numerous strain-specific genes, transposases, IS elements, and genes previously identified as being differentially expressed under aggregation conditions (11). In two genomes (F. nucleatum subsp. nucleatum ATCC 25586 and F. nucleatum subsp. polymorphum ATCC 10953), a synteny break (indicated with a black line) placed the canonical fadA gene (FN0264) in this region rather than a second radD gene. Several regions similar to this are located in the genome.

To gain evidence for whether the active invader-specific orthogroups were acquired in the adaptive radiation (Fig. 1), we examined the occurrence of the 44 active invader-specific orthogroups in the closely related outgroup Leptotrichia buccalis, the presence of which would suggest that the Fusobacterium ancestor possessed them prior to the radiation and that they were lost by passive invaders. A minority (17 of 44 [39%]) of active invader-specific orthogroups were found in L. buccalis. Those predicted from this analysis to have been lost in the radiation of passive invaders included genes encoding a membrane-bound two-component sensor/histidine kinase, a META domain-containing protein, a hypothetical transporter, and a secreted peptidoglycan catabolism protein. Genes not present in the outgroup and likely to have been acquired by active invaders in the radiation included four of the six orthogroups previously implicated in virulence in other species, including those encoding branched-chain amino acid transport and type IV pili (40, 41, 44–46), as well as 15 of the 19 genes encoding predicted membrane-bound or secreted proteins, highlighting acquired changes in surface organization of the active invaders.

Interestingly, the complement of MORN2 domains remained relatively unchanged between Leptotrichia and the passive invaders but was greatly expanded in the active invader clades (Fig. 2D). In contrast, FadA and BMC domain-containing genes, both present in Leptotrichia, appeared to have been lost in the passive invaders, as well as in F. mortiferum. YadA was not observed in Leptotrichia but may not be necessary for active invasion, since the active invader clade C (F. ulcerans/F. varium) lacked genes containing these domains. That F. mortiferum contained no YadA, FadA, or BMC domains and only a modest expansion of MORN2 domains suggested that it might not be capable of active invasion. Overall, our data predict an evolutionary model whereby (i) the last common ancestor of Leptotrichia and Fusobacterium shared active invasion-associated FadA, BMC, RadD, and MORN2 domain-containing genes, (ii) FadA and BMC were lost during the evolution of F. mortiferum, F. necrophorum, and F. gonidiaformans, and (iii) MORN2- and, to some degree, RadD domain-containing genes underwent expansion in the active invaders. In support of this model, we observed that at least three Leptotrichiaceae genomes in the Pfam database also shared FadA, RadD, and BMC domains and a small set of ancestral MORN2 domains (data not shown).

MORN2 evolution and function.Of the expanded Pfam domains, MORN2 was among the most intriguing because it represented the most frequent domain in active invader genomes (115 to 250 MORN2 domain copies per genome), and nothing was known about the function of proteins containing this domain. The MORN2 domain itself is 22 to 23 amino acids long, often found in multiple copies per gene, and was highly variable across the strains in our study. MORN2-containing proteins are incredibly diverse in the organization and grouping of domains, which is comprehensively illustrated at http://pfam.xfam.org/family/PF07661. The number of genes containing MORN2 domains (MORN2 genes) was also variable, with members of the same subspecies differing in their MORN2 gene content by as much as 25% (Fig. 4). The small number of MORN2 genes found in the passive invaders was also present in nearly all active invaders, suggesting that MORN2 genes evolved and expanded in active invader species from an ancestral set of MORN2 genes held by the last common ancestor at the time of the adaptive radiation. Supporting this, we observed (i) no evidence that MORN2 genes were recently acquired on genomic islands (using Islandviewer [56]), (ii) MORN2 genes clustered in active invader genomes (Fig. 3), with evidence of recent local duplication in some species (Fig. 5), and (iii) there was significantly greater sequence identity between pairs of MORN2 domains that were spatially close but not in the same protein (38.3%), compared to those >10 kb away (31.3%) (P < 2.2e−16). We also observed an unusually high level of overall chromosomal rearrangements in actively invading Fusobacterium strains, which could be linked to the high rate of genome rearrangements within MORN2 regions (see Fig. S3 and Text S1 in the supplemental material).

• Open in new tab
• Download powerpoint

FIG 4 

Species distribution of MORN2 orthogroups. Each column represents an orthogroup. The boxes are colored according to the average number of MORN2 domains per protein, and for each orthogroup containing two or more paralogs, the numbers indicate the number of paralogs. There is tremendous variation in the number and structure of proteins containing MORN2 domains across genomes. Active invaders (species names shaded in dark gray) contain longer, more complex MORN2-containing proteins than passive invaders (species names shaded in light gray). Passive invaders contain only a small complement of short, “ancestral” orthogroups, which are present across all species. Expansions of active invader orthogroups, especially in F. periodonticum, can be seen by the presence of paralogs.

• Open in new tab
• Download powerpoint

FIG 5 

Expansion of an “ancestral” MORN2 orthogroup in F. periodonticum. This orthogroup is present in all species, including the passive invaders, but it is greatly expanded in F. periodonticum. In F. periodonticum, this expanded cluster is also located near other known virulence-related adhesins, including RadD and YadA family members. All members of this orthogroup are present in this small region of the genome. Arrows represent syntenic orthogroup members. The orthogroup identification no. for this expanded group is 844916212 (see Table S2 in the supplemental material).

Since nothing was previously known about the function of MORN2 genes, we sought clues to their function by comparing the ecologies and lifestyles of other bacteria that exhibited expanded sets of MORN2 domains. Remarkably, although MORN2 domains were present in a diverse set of bacterial genomes (see Text S1 in the supplemental material), the only other bacterium that contained >100 MORN2 domains, like active invader Fusobacterium, was Helicobacter bilis ATCC 43879, which colonizes the bile, liver, and intestines (24). The closely related species Helicobacter hepaticus also possessed a relatively large number of MORN2 domains, 28 copies, which was higher than the number in the passively invading fusobacteria. H. bilis and the related H. hepaticus are both involved in colitis, hepatitis, and, similar to actively invading Fusobacterium spp. exhibit carcinogenic potential and association with preterm birth in animals (24–26). Although it is unknown whether H. bilis or H. hepaticus actively invades host cells, given (i) the striking similarity in pathology associated with these two species and actively invading Fusobacterium and (ii) the fact that MORN2 genes colocalize with the expanded FadA family and RadD family genes known to be involved in adhesion and virulence in Fusobacterium, we propose a testable model that predicts a related role for MORN2-containing proteins in adhesion and pathogenicity and possibly in promoting cancer or cancer progression.