Clostridioides difficile is a major worldwide cause of antibiotic-associated gastrointestinal infection. Two toxins (TcdA and TcdB) are responsible for C. diffile pathogenicity, but genetic variants within these toxins complicate the development of broad-spectrum diagnostics, therapeutics and vaccines. Here we provide a global classification and analysis of available C. difficile toxin sequences and introduce a new open online database (diffbase.uwaterloo.ca) to serve the unmet needs of the clinical and research community. Our analysis partitions TcdA and TcdB genes into 7 and 12 distinct groups which provides a new method for sequence-based C. difficile toxin subtyping. Our analysis revealed that recombination has driven extensive diversification of TcdB in particular, resulting in TcdB subtypes with distinct antigenic, functional, and phenotypic properties. As validation of our method, we were able to rapidly subtype a new dataset of 351 clinical strains from Brigham and Women’s Hospital, predicting their phenotypic and clinical features. Lastly, based on sequence analysis, we identified conserved regions in TcdB that represent ideal targets for the development of universal C. difficile therapeutics and diagnostics.
Clostridioides difficile (formerly Clostridium difficile ) is a diverse group of Gram-positive spore-forming anaerobic bacteria [1]. Toxigenic strains have become important opportunistic pathogens to humans. Their spores are widespread and can colonize human and animal colons after disruption of the gut microflora, most notably due to antibiotic treatment. C. difficile infection (CDI) results in a range of symptoms from self-limiting diarrhea to severe pseudomembranous enterocolitis and death [2–7]. It is the most frequent cause of healthcare-associated gastrointestinal infections across developed countries worldwide [2–5,8].
Ribotyping (RT), which compares intergenic spacers between ribosomal RNA genes, is widely utilized to categorize C. difficile linages [5,9]. Various other methods including multilocus sequence typing based on allelic variation of housekeeping genes and whole genome sequencing analysis have also been adopted to further discriminate strains [5,9–13]. Phylogenetic analyses revealed a growing diverse population [1,14–16]. In recent years, there is an emergence and spreading of various epidemic hypervirulent strains such as the RT027 clonal lineage, which first caused outbreaks in 2000–2003 in North America and is associated with increased disease severity and mortality [17–20]. RT078 is an emerging hypervirulent lineage which is also the dominant type found in domesticated animals [21,22]. There are also geographic differences, for instance, RT017 has become a dominant lineage in Japan and Korea [23].
The major virulence factors in toxigenic C. difficile strains are two homologous large protein toxins, TcdA (~300 kDa) and TcdB (~270 kDa) [24–27]. Nontoxigenic C. difficile strains without these toxins exist and can colonize humans and animals, but do not cause diseases [28]. TcdA and TcdB share overall ~66% sequence similarity and belong to the large clostridial toxin (LCT) family, which include TcsH and TcsL in Paeniclostridium sordellii, Tcnα in Clostridium novyi, and TpeL in Clostridium perfringens [5,6,8,9,24,25,29–31]. TcsH and TcsL can be considered orthologs of TcdA and TcdB, respectively, with TcsH sharing ~77% sequence identity with TcdA and TcsL sharing ~76% identity with TcdB [32] (S1 Fig). TcdA and TcdB share a protein domain architecture consisting of an N-terminal glucosyltransferase domain (GTD), followed by a cysteine protease domain (CPD), an intermingled membrane translocation delivery domain and receptor-binding domain (DRBD), and a large C-terminal combined repetitive oligopeptides domain (CROPs) (S1 Fig). After binding, endocytosis, and translocation across endosomal membranes into the cytosol of host cells, these toxins glucosylate and inactivate host Ras/Rho family of small GTPases, leading to disruption of the actin cytoskeleton, cell rounding, and ultimately cell death [33].
TcdA and TcdB were first identified, characterized and sequenced in the 1980s [34–38] and the toxin sequences from a reference strain (VPI10463) have been widely used as the standard in diagnostic and therapeutic development. However, sequence variations in the toxin genes exist across C. difficile strains and could affect receptor-binding specificity, preferences toward distinct small GTPases, overall toxicity, and antigenicity. For instance, strains such as R20291 (belonging to RT027) produce a TcdB variant with ~8% of residue differences from the reference TcdB, which exhibited a significant impact on its immunogenicity: mice immunized with the reference TcdB developed resistance to the same TcdB, but all died when challenged with this variant TcdB [39], and several antibodies raised against the reference TcdB, including the FDA approved therapeutic antibody bezlotoxumab, either do not recognize or have lower efficacy against this TcdB variant [39–41]. Furthermore, this TcdB variant also lacks the ability to recognize frizzled (FZD) proteins, which are one of the major receptors for the reference TcdB, due to residue changes at the FZD-binding interface [40,42–45].
These toxin variations pose a significant challenge for developing effective broad-spectrum diagnostic assays, therapeutic antibodies, and vaccines. Understanding variations in toxins is a key step to address this challenge and may also reveal their potential evolutionary paths and functional differences. A toxinotyping method has been previously developed utilizing PCR-based amplification of toxin gene fragments and analyzing polymorphism with restriction enzyme digestions, which can distinguish over 34 toxinotypes [46,47]. Although toxinotyping highlights the variation among toxin genes, it lacks the resolution to understand the molecular basis for diversification of toxins and sequence-function relationships.
Rapid growth of genomic sequencing of C. difficile strains in recent years provides an opportunity to analyze and categorize the diversification of TcdA and TcdB with single residue resolution. Here we performed a comprehensive analysis of nearly all available C. difficile TcdA and TcdB sequences, including assembly and analysis of 6,492 new genomes, with the goal to 1) build a comprehensive central database of C. difficile toxin sequences; 2) better understand the mechanisms underlying TcdA and TcdB diversification; and 3) develop a system to classify TcdA and TcdB into subtypes that allow clinicians and researchers to categorize and predict functional-immunological variations of any future sequenced C. difficile isolates.
To build a comprehensive database of TcdA and TcdB sequences, we combined data from NCBI GenBank and the NCBI short-read archive (SRA). From 2,347 C. difficile genomes in GenBank, we identified an initial set of 1,633 tcdA and 2,028 tcdB genes. We then developed a computational pipeline for automated retrieval of C. difficile genomes from the SRA, de novo genome assembly, genome annotation, and extraction of tcdA and tcdB genes (see Methods). Using this pipeline, we assembled the genomes of 6,492 C. difficile isolates and identified an additional 3,542 tcdA and 5,994 tcdB genes (Table 1). Combining both sources, we identified 5,175 TcdA and 8,022 TcdB encoding sequences.
We then carried out alignments of all toxin protein sequences. The TcdB alignment covered the entire sequence (1–2366), with 712 (30%) of the positions showing variations across all domains. The TcdA alignment possessed much lower variation than TcdB within the 1–1874 region as it had only 168 (9%) variable sites, but its CROPs domain (1831–2710) contained an extremely high degree of variation in the number and length of repeats: from 3 repeats in the shortest variant to 45 in the longest variant, and 32 in the reference TcdA variant from VPI 10463 (S2 Fig). This is likely generated by homologous recombination due to the repetitive nature of this region. The CROPs domain is composed of long repeats (LRs) of ~30 residues and short repeats (SRs) of ~19–24 residues [27]. The CROPs domain in TcdA is not only repetitive at a protein sequence level, but also showed a high degree of repetitiveness at a DNA level, whereas the repetitiveness of the CROPs domain in TcdB is largely limited to the protein level [48,49], which may account for frequent recombination in TcdA-CROPs but not in TcdB-CROPs. It is also important to note that the repetitive nature of the TcdA-CROPS region at the DNA level may result in assembly errors, which may inflate the apparent variation in this region.
In total, there were 116 unique TcdA protein sequences and 212 unique TcdB protein sequences. We then clustered these sequences into distinct subfamilies (“subtypes”) using average linkage hierarchical clustering (see Methods). Analysis of TcdB is based on full-length sequences, but TcdA is limited to the 1–1874 region to avoid the highly variable CROPs domain. In addition, we also included TcsH and TcsL sequences in our analysis. Clustering produced 7 distinct TcdA subtypes which we labeled A1-A7, and 12 distinct TcdB subtypes which we labeled B1-B12 (Fig 1). B1-B4 were ordered to match previous literature [50–52], and additional subtypes were ranked based on total frequency of occurrence in GenBank and NCBI-SRA. Each unique sequence was then further numbered following a period within its subtype (e.g. B1.1, B1.2, B1.3, etc.). Sequences within the same TcdA and TcdB subtype demonstrate strong pairwise similarities, and weak similarities between subtypes (Fig 1A and 1D). Quantitative analysis revealed that thresholds of 99.4% (TcdA) and 97% (TcdB) can be used to effectively assign toxin sequences to these subtypes (S3 Fig). We then selected one representative sequence for each subtype and carried out phylogenetic analysis and pairwise comparison. TcdA subtypes A1 to A6 possess higher similarities (>97.9%) and clustered together, with A7 forming a divergent lineage (Fig 1B and 1C). A7 is a unique sequence with only 85.3% to 85.6% identity to others (Fig 1C). The entire TcdA family was further outgrouped by TcsH as expected (Fig 1B). TcdB also formed a monophyletic family that was outgrouped by TcsL and a second lineage of TcsL-related proteins (Fig 1E). TcdB subtypes can be subdivided into three groups, one including B3, B4, B7, B8, a second including B2, B9, B10, B11, and a third including B1, B5, B6, and B12 (Fig 1E). The lowest identity among TcdB subtypes is 85.3% (between B4 and B12, Fig 1F). A7, B10, B11, and B12 represent rare divergent subtypes recently reported: A7 is in strain RA09-70, which does not express TcdB [53]; B10, B11, and B12 were identified recently from strains CD10-165, CD160, and 173070, respectively [53], and all three strains lack TcdA.
By mining unassembled C. difficile genomes from the SRA, we were able to discover 125 TcdA and TcdB protein sequences that were not represented in GenBank. Most novel toxin variants clustered into subtypes A1 (N = 25), B1 (N = 52), and subtypes A2 (N = 10) and B2 (N = 12) (S1 Table). However, three highly divergent TcdA variants identified from SRA datasets formed new subtypes not represented in GenBank. These include subtypes A4 from strain L;13.7548369.T (SRS1486661), A5 from strain ECDC-088 (SRS1486236), and A6 from strain ECDC-009 (SRS1486256), all of which are clinical isolates. All three of these strains contained truncated/partial TcdB variants which represent putative pseudogenes.
To link our subtyping with known clinical C. difficile strains, we manually curated subtype assignments for a set of 70 C. difficile strains selected from the literature, which covers known toxinotypes, and compared subtypes with toxinotypes, ribotypes, and whether the strain produces the third toxin known as C. difficile transferase toxin (CDT) (S2 Table). The majority express an A1/B1 subtype combination and include reference strains 630 and VPI 10463 that express the widely used standard TcdA and TcdB sequence (defined as A1.1 and B1.1, S2 Table). The group that expresses a combination of A2/B2 is the second largest and includes hypervirulent RT027 strains R12087 and R20291. The group expressing A3/B5 include strains (e.g. M120 and NAP07) classified as RT078. Subtype B3 is mainly expressed in strains (e.g. 1470) belonging to RT017, which lacks TcdA. Other pairings in the table include A2/B9, A3/B6, A2/B7, and A1/B3. The table includes many strains that do not express functional TcdA, which can express B1, B2, B3, B4, B5, B7, B8, B10, B11, or B12; one strain that only expresses TcdA but not TcdB (A7 in RA09-70); four strains that only express CDT; and one strain (SLO037) that expresses none of the three toxins. This table represents a small portion of C. difficile strains and a full list of a total 1,640 C. difficile strains from the NCBI database with their toxin subtypes noted is included as S3 Table.
In general, phylogenetic subtyping of C. difficile toxins correlated well with previously identified toxinotypes, but at greater resolution by analyzing TcdA and TcdB separately (S2 Table, see Discussion). There was less congruence with ribotypes, however, as different subtypes were found in the same ribotype strains, and the same subtype was found in different ribotype strains. Therefore, neither toxinotype nor ribotype were able to accurately categorize toxins based on phylogenetic relationships (S2 Table). Subtyping was capable of capturing the full phylogenetic diversity of TcdA and TcdB available in previously known and new strains.
To evaluate the phylogenomic distribution of toxin subtypes across C. difficile, we constructed a whole-genome based phylogeny of 1,934 complete C. difficile genomes based on 14,194 SNP positions across 88 conserved marker genes (Fig 2A and S3 Table) (see Methods). The genome tree is highly consistent with known phylogenetic relationships, as the previously identified clades 1–5 are represented by distinct lineages [14] (Fig 2A). Two of the three divergent environmental lineages C-I and C-II are also present as divergent branches (Fig 2A).
A total of 1,640 (84.8%) C. difficile strains were found to encode TcdA and/or TcdB, while the remainder (294, 15.2%) lack toxin genes. The predicted toxin subtypes across the C. difficile genome tree demonstrate strong clade associations, and therefore are highly congruent with strain phylogenetic relationships. The congruency between subtype and phylogeny provides further support for our toxin classification (Fig 2A). For example, subtype A1/B1 which includes reference strains 630 and VPI 10463 is most common among toxin-containing strains (979, 59.7%) and associated with clade 1 (Fig 2B). A2/B2 was second most common and associated with clade 2, A3/B6 with clade 3, -/B3 with clade 4, and A3/B5 with clade 5. Also prevalent were types -/B1, A1/-, and A2/B7 (Fig 2B). Deviations from the A1/B1 toxin type are often associated with the emergence of numerous hypervirulent and epidemic outbreak strains such as A2/B2 (RT027), A3/B5 (RT078), and -/B3 (RT017) (Fig 2A). Interestingly, the highly divergent environmental lineages encode the highly divergent TcdB subtypes B10 and B11 (C-I) and B12 (C-II) (Fig 2A and S2 Table). This is consistent with an early divergence of B10-B12 in C. difficile evolution, predating the emergence of TcdB subtypes found in the other clinical strains.
Interestingly, we also observed rare lateral transfer events involving only one of the two toxin genes to create hybrid strains containing new subtype combinations. Examples include the spontaneous emergence of an A1/B3 strain within clade 1, and the emergence of an A1 /B2 strain in clade 2 (Fig 2A). Thus, through lateral transfer and homologous recombination, subtype B3 has likely replaced B1 in a clade 1 strain, and subtype A1 has likely replaced A2 in a clade 2 strain. Furthermore, we observed many independent clades containing tcdA-/tcdB- C. difficile strains (e.g., see six lineages marked by asterisks in Fig 2A). This is consistent with previously reported “defective” toxin clades [54], and indicates numerous independent losses of the pathogenicity locus throughout C. difficile evolution.
As an independent test dataset for our toxin subtyping method, we examined 351 genomes of C. difficile isolates derived from a clinical cohort from Brigham and Women's Hospital (BWH) in Boston (Fig 2C) [55]. As they were not included in our initial database, they are ideal for testing the robustness and effectiveness of our subtype classification. All identified toxins could be accurately assigned to our reference sequences, with most (97%) aligning with 100% identity to our database, and the remainder aligning with > = 99.8% identity. Out of 351 total strains, 62 (17.7%) were toxin deficient, while 289 (82.3%) contained TcdA and/or TcdB genes (S4 Table). Of these, there were 12 distinct subtype combinations, with frequencies similar to those observed in the NCBI dataset. A1/B1 strains were most common (N = 222), followed by A2/B2 (N = 24), -/B3 (N = 11), and A3/B5 (N = 10) (Fig 2C). Therefore, our method was able to rapidly and automatically classify a large dataset of 351 clinically relevant C. difficile isolates, with all sequences represented in our current classification.
We next focused on understanding the evolution of TcdA and TcdB variants and mechanisms for their diversification. To visualize global patterns of variation within TcdA and TcdB, we developed a haplotype coloring algorithm (https://github.com/doxeylab/haploColor) based on previous methods for genome visualization [56]. First, sequences are painted black where they matched the reference sequence (i.e., B1.1). Then, remaining positions were painted different colors where they matched selected other subtypes (Fig 3A): blue when matching B5.1, gold when matching B7.1, and green when matching TcsL. The result of this algorithm applied to the TcdB alignment revealed a striking block-like and highly mosaic pattern of amino acid variation, which strongly indicates recombination between subtypes (Fig 3A). B1, B5, and B6 are composed of a B1-like variation (black) pattern across their full-length sequences, while B4 and B7 are composed of a B7-like pattern (gold) across their full-length sequences. B2, B3, B8, and B9, however, possess a mosaic combination of B1-like and B7-like patterns. B3, B4, B7, and B8 share a distinct B7-like pattern of amino acid variation across their N-terminal region including the GTD and CPD domains, but when examining the DRBD, the B7-like pattern is shared by a different set of subtypes (B2, B4, B7, B10, and B11). These patterns indicate ancestral within-gene (“intragenic”) recombination events involving distinct regions of TcdB. As a statistical test of recombination, we further performed phylogenetic network analysis using SplitsTree [57]. Consistent with patterns of amino acid variation and per-domain phylogenetic analysis, network analysis revealed significant evidence of recombination within TcdB (p = 0; Phi test for recombination) (S4 Fig). In contrast to TcdB, TcdA (1–1874) produced homogeneous patterns of variation across each subtype (S5 Fig) and did not display evidence of recombination in network analysis (p = 0.186) (S4 Fig), indicating that recombination occurs frequently only in TcdB, but not in TcdA.
We further performed separate phylogenetic analyses of each domain (GTD, CPD, DRBD, and CROP) of TcdB (Fig 3B). The phylogenetic tree of each domain produced two main groups (labeled i and ii), which correspond with the B1-like and B7-like patterns revealed in the alignment visualization (Fig 3A). Each subtype can therefore be described as a chimeric combination of type “i” (B1-like) or type “ii” (B7-like) domains (Fig 3C). Based on the per-domain phylogenetic relationships and recombination patterns, we formulated a potential evolutionary model for the origin of TcdB subtypes (Fig 3D). An early TcdB ancestor split into two main groups: (i) B1, B5, and B6; and (ii) B4 and B7. Subtype B2 likely originated by a recombination event fusing an ancestral type i and type ii toxin. B9 likely originated from a recombination event between B1 and B2, B3 from a recombination event between B1 and a type ii toxin, and B8 from a recombination event between B5 and a type ii toxin. Subtypes B10-B12, which are recently identified rare variants, are early diverging lineages since they consistently outgrouped other subtypes in phylogenetic analysis (Fig 3B), consistent with their divergent lineages among other strains (Fig 2A).
In addition to these major ancestral recombination events, we also identified a considerable degree of “microrecombination” events involving exchange of small segments between subtypes. For example, a single TcdB sequence (B1.59) from subtype B1 has acquired an N-terminal segment that is clearly derived from subtype B2 or B9 (Figs 3A and S6). This unique TcdB gene, which appears to be the result of a spontaneous recombination event between a B1 and B2-containing strain, is derived from a newly assembled clinical isolate from a Fidaxomicin clinical trial (SRS1378602). A second similar example is B1.58 from a clinical isolate (ECDC-040, SRS1486176), which has acquired a DRBD and CROPS segment from a B2-containing strain (S6 Fig). Fourteen such cases of microrecombination including these are depicted in S6 Fig. TcdB in particular appears to have diversified through an extensive degree of intragenic recombination involving both large and small segments.
The value of subtyping classification is to facilitate a molecular understanding of the impact of sequence variations on function and antigenicity. For instance, our sequence alignment divides the GTD into two groups: one contains B3, B4, B7, and B8; and the rest form another group (Figs 3B and S7). Previous studies have reported two types of cell-rounding effects: TcdB1 and B2 are known to induce rounded cells with many protrusions remaining attached to cell culture plate, whereas TcdB from the strain 1470 and 8864 have been reported to cause rounded cells without protrusions, which is similar to TcsL [58]. It has been proposed that this is a result of the altered specificity of their GTD in targeting different small GTPases [32,58]. TcdB in strain 1470 is classified as B3, and the strain 8864 expresses B4, thus our classification predicts that the group containing B3/4/7/8 induces a TcsL-like cell rounding phenotype. This is indeed the case for two recently reported clinical strains HSJD-312 and HMX152: both express toxins classified as B7 under our subtyping system (S2 Table) and have been reported to induce TcsL-like cell rounding [59].
Another well-characterized functional motif in TcdB is its FZD-binding interface, with key residues clearly defined by the co-crystal structure [42,43,60,61]. It has been reported that B2 lacks the ability to bind FZDs due to residue variations at FZD-binding interfaces [40,44,45]. To survey whether these variations may also exist in other subtypes, we aligned the key residues across all TcdB sequences and visualized them. As shown in Fig 4A, the FZD-binding motif is highly conserved across B1/3/5/6/8/9, while B2/4/7/10 share the same set of residue changes. Thus, B4/7/10 are predicted to lack FZD-binding capability similar to B2. B11 contains a subset of residue changes found in B2 within this region and likely also has reduced binding to FZDs. This pattern is consistent with the phylogenetic alignment of the DRBD domain, in which B1/3/5/6/8/9 form group i and B2/4/7/10/11 form group ii (Fig 3B and 3C). Interestingly, although most B2 variants possess FZD-binding site substitutions, there are a few exceptions that contain a largely intact FZD binding site. In particular, B2.12 assembled from strain 2007223 (ERS001491) contains only a single amino acid substitution (F1597S) in this region. Examination of the alignment reveals that this is likely due to a microrecombination event that has replaced most of the FZD binding site with a B1-like segment (Figs 4A and S6). A similar scenario occurred in a member of subtype B7, in which a B1-like segment has partially replaced this region (S6 Fig).
Sequence variations between subtypes could also have a drastic impact on the efficacy of therapeutic antibodies and vaccines. Bezlotoxumab from Merck is the only monoclonal antibody against TcdB that has so far been approved by the FDA and is currently used to reduce the recurrence of CDI [62]. This antibody was generated using fragments of TcdB1 as antigens and its epitope sites (located at the N-terminal of CROPs) have been established through crystallography [63]. We thus aligned all key residues within its epitope across all TcdB sequences, which revealed extensive residue changes largely conserved in B2/4/7/9/10/11 (Fig 4A). This is consistent with our alignment of the CROPs domain that group B2/4/7/9/10/11 together (Fig 3B and 3C). It has been shown that bezlotoxumab exhibited as low as ~200-fold reduction in neutralization efficacy against TcdB from several RT027 strains, which likely express B2, compared with its efficacy against B1 from VPI10463 [41,64]. It also showed a similarly low efficacy against a strain 8864, which expresses B4. These results indicate that bezlotoxumab does not have good efficacy against CDI caused by strains that express B2/4/7/9. Furthermore, there are also a few amino acid changes within the epitope region in B5/B8, and it has been shown that bezlotoxumab has ~60-fold reduction in efficacy against the TcdB from a RT078 strain [41], which likely express B5 (S1 Table). These results clearly indicate that subtype classification of toxins will be able to guide the use of bezlotoxumab in the clinic.
In addition to bezlotoxumab, we also examined another monoclonal antibody PA41, which is under development [41], and a single-domain antibody (also known as VHH or nanobody) E3 [65]. The epitopes for both have been well established through co-crystal structures [65,66]. Both recognize the GTD domain, with E3 recognizing the N-terminus of TcdB (Fig 4A). The epitope site for PA41 is highly conserved across most subtypes except a single residue change (Y323H) in B3. This is consistent with the previous finding that PA41 can potently neutralize TcdB from many different strains except RT017 strains, which express B3 [41]. The epitope site for E3 is conserved in most subtypes except a single residue change (I58T or A) in B3/4/7/10/11/12, and the impact of this single residue change remains to be examined experimentally.
Finally, we mapped evolutionary conservation across all available TcdB sequences onto the recently reported crystal structure of TcdB [65] (Fig 4B). Relative to the reference TcdB1 sequence, amino acid variants are common across the full-length TcdB sequence and occur throughout each domain (Fig 4A) but some regions (e.g., N-terminus of the GTD, C-terminus of CROPS domain, segments of the pore-forming region of the DRBD and C-terminus of the CROPS domain) were highly conserved. Based on structure, we identified eight conserved surface patches containing universally conserved residues which represent potential key therapeutic targets for developing broad-spectrum diagnostics, antibodies, and vaccines (Fig 4B).
To address the needs of the research and clinical community in understanding toxin subtyping and variations, we developed an online open database freely accessible at diffbase.uwaterloo.ca . DiffBase stores all unique TcdA and TcdB sequences identified to date from the NCBI and SRA and organizes sequences into our subtype classification scheme. Different subtypes and individual sequences can be explored and visualized in reference trees, with additional information such as source strains, and links to other resources (S8 Fig). In addition, users can query their own TcdA or TcdB sequences against the database using a built-in BLAST interface, which will report the top matching sequences in the database and provide toxin classifications and other related information. To keep up with new sequences and information concerning TcdA and TcdB, DiffBase facilitates community feedback and allows users to submit new information to be added to the next iteration of the database.
Here we created the largest database to date capturing available TcdA and TcdB sequence diversity. This up-to-date collection includes genes from sequenced C. difficile isolates in GenBank, as well as thousands of genomes that were assembled, annotated and analyzed from the NCBI short-read archive. We clustered TcdA and TcdB variants into phylogenetic subtypes, which provided a robust classification that is both congruent with C. difficile genome phylogeny as well as variation in functional and therapeutically relevant amino acids including TcdB regions targeted by existing monoclonal antibodies. Our analysis revealed that TcdB undergoes extensive homologous recombination, and its potential evolutionary history is proposed based on recombination among various subtypes. Finally, our analysis revealed eight conserved patches across the TcdB structure, which will facilitate future studies that aim to develop “universal” C. difficile therapeutics that broadly target all TcdB subtypes.
In general, there is some agreement between previously defined toxinotypes and our toxin subtypes, but subtyping provides additional information as it is able to describe TcdA and TcdB separately. For example, toxinotype 0 associates largely with the A1/B1 subtype, toxinotype III associates largely with A2/B2, toxinotype IV with A3/B6, toxinotype VIII associates with -/B3, toxinotype IX associates largely with A2/B7, toxinotype X associates with -/B4, and so on (S1 Table). However, subtype A3/B5 associated with toxinotypes V, VI, VII, XVI, XXVIII, all of which are found in clade 5 strains. Moving forward, with improved abilities to perform genome sequencing of clinical isolates, it will be increasingly possible to classify strains based on their genome-wide phylogenetic relationships as well as their toxin subtypes.
In comparison to TcdA, our analysis identified a much greater degree of sequence variation within TcdB and a larger number of subtypes. Given that we see evidence for extreme recombination in TcdB but not TcdA, it is possible that there is a greater selective pressure for positive selection and diversification of TcdB. We speculate that intragenic recombination of TcdB may drive antigenic diversification, whereas in TcdA this process may be driven by truncation and variation of its C-terminal CROPS region. The CROPs domain showed similarity with carbohydrate-binding proteins and may contribute to toxin attachment to cells by binding to carbohydrate moieties (27, 42, 44). The CROPs domain may also act as a chaperone that protects other domains (45). Possibly due to its repetitive nature, the CROPs domain is often the region that induces strong immune responses. It remains to be determined whether frequent recombinations/changes in TcdA-CROPs may alter its function and/or antigenicity.
These findings further suggest that TcdB may play a central role in C. difficile pathogenesis, which is consistent with previous findings that TcdA-/TcdB+ mutant C. difficile strains are fully virulent, whereas TcdA+/TcdB- strains are attenuated in multiple mouse models [26,67]. It has also been suggested that TcdB is the primary factor for inducing the host immune and inflammatory responses in mouse models [26]. The key role of TcdB in CDI is further confirmed by the findings that an antibody that neutralizes TcdB (bezlotoxumab), but not another one that neutralizes TcdA (actoxumab, Merck), conferred protection against CDI in gnotobiotic piglets [68] and reduced CDI recurrence in humans [62,69] and it is also consistent with the fact that many clinical isolates only express TcdB [70]. An exception to a dominant role for TcdB is the very rare TcdA+ TcdB- strain. It is noteworthy that one such strain identified in GenBank contains the single most divergent TcdA sequence (subtype A7) [53], which may have diverged to acquire a pathogenic functionality without requiring TcdB.
For such recombination events to have occurred in TcdB, it is possible that phylogenetically distinct C. difficile strains containing different toxin subtypes coexisted within the same host individuals, exchanged genetic material and recombined to produce new recombinant forms. Co-infection with different C. difficile ribotypes has been recently reported in a clinical case study [71]. Theoretically, co-infection does not need to occur frequently to promote recombination. A single individual containing two or more C. difficile strains may be sufficient to promote recombination, generating hybrid toxins with different regions derived from different subtypes. The new recombinant strain can then increase in frequency through genetic drift and/or selection, as well as through transmission to other individuals. Our analysis suggests that this process has not only occurred frequently in the past as a mechanism by which different subtypes originated, but that it may be a frequent and ongoing process in new clinical isolates (e.g., B1.59 from SRS1378602). Consideration of intragenic recombination and how it may shape TcdB function and toxicity will be important in efforts to understand the emergence of new C. difficile hypervirulent strains and develop targeted therapeutic interventions.
Recombination offers considerable adaptive benefits to proteins by facilitating rapid mutation of a sequence by exchange of entire segments as opposed to the relatively slower process of single point mutations. In this way, proteins can diversify by shuffling a few basic building blocks such as protein domains. In pathogens, recombination plays a major role in pathoadaptive evolution by facilitating rapid “switching” of virulence factors and antigenic proteins [72,73]. Antigenic recombination can promote the sudden avoidance of immune recognition (antigenic escape), which has been demonstrated for the C. difficile S-layer gene [74]. In the case of TcdB, intragenic recombination may generate new hybrid toxins composed of different domains types and functions. In theory, recombination could also generate resistance to therapeutics by replacing entire binding interfaces with compatible regions from other toxins that possess drug-resistant mutations. Recombination-mediated domain shuffling not only describes TcdB sequence patterns and phylogenetic relationships, but also provides an explanation for important functional differences between TcdB variants. For example, the exchange of a B7-like GTD between subtypes B3, B4, B7 and B8, correlates with the TcsL-like clumping and rounding phenotype. Also intriguing are the many microrecombination events that have occurred in the DRBD region which overlap with FZD-binding site. For example, likely due to partial homologous recombination with a B1-like toxin, one B2 variant (B2.12, strain 2007223 from ERS001491) contains an intact FZD-binding interface with only a single amino acid substitution (F1597S). This suggests that intragenic recombination in TcdB may promote rapid evolutionary switching between receptor-binding activities or affinities.
Given the extent of TcdB diversification and its primary role in virulence, it is critically important to identify conserved regions that can be targeted for therapeutic and diagnostic applications. Sequence conservation mapped to protein structure also revealed at least 8 distinct surface patches containing a high density of universally conserved residues across all TcdB subtypes, which represent promising regions for the development of inhibitors. Importantly, the binding site for the antibody therapeutic bezlotoxumab, which is commonly used to treat C. difficile infections, was not among these and instead displayed considerable variation across TcdB subtypes with B2, B4, B7, B9, B10, B11, in particular displaying 7–8 likely destabilizing substitutions. Although the common B1 subtype of TcdB is largely conserved across this region, based on analysis, it is possible that intragenic recombination with other strains (e.g., a B2-containing strain) could generate spontaneous resistance to bezlotoxumab by replacing this region with a B2-like segment. Future efforts to target highly conserved clusters of surface-exposed residues on the TcdB structure may yield promising candidates for therapeutic or vaccine development.
Finally, based on sequence-based classification of tcdA and tcdB genes, we propose a revised scheme for naming these genes in future studies. In this scheme, a newly identified TcdA or TcdB sequence may be aligned to our reference database and named based on the top hit according to sequence identity. In order to enable automated subtyping of new tcdA and tcdB genes and facilitate community collaboration and data sharing, we have developed a freely available online database (DiffBase) for use by the C. difficile clinical and research community. In the future, clinicians will be able to query toxin sequences from clinical isolates and immediately determine the toxin subtype, which will help them decide on therapeutic strategies. For instance, among the 351 clinical cases in the BWH dataset, there are 34 cases expressing B2/B4/B7/B9, for which treatment of bezlotoxumab would not be effective. Therefore, toxin subtyping will guide proper choices of clinical treatment in consideration of toxin variations and allow researchers to monitor the ongoing evolution and diversification of C. difficile.
Note: Related work was published during the preparation of this manuscript that performed a subtype analysis of TcdB [50]. The analysis was limited to TcdB only and also included a smaller set (3,269) of genomes, resulting in fewer (8) TcdB subtypes. Their classification is consistent with our analysis, but does not capture the full diversity of available sequences or subdivide sequences according to a comprehensive analysis of recombination patterns. Furthermore, no online database was provided to facilitate automated subtyping analysis. Our work is also unique through its model of toxin evolution, the analysis of conserved and variable regions and their impact on antigenicity, the discovery that the current therapeutic antibody is only effective on a subset of toxins, and its analysis of a validated clinical dataset from Brigham and Women’s Hospital.
A set of Clostridiodes difficile sequencing runs was retrieved from the NCBI short read archive (SRA) by text query for “Clostridioides difficile ” on June 20th, 2019. Metagenomic samples were omitted, leaving only genomic samples to reduce the chance of contamination from other bacterial species. Sequencing runs were downloaded using the fasterq-dump module of the SRA toolkit. To account for multiple library preparation methods and adapters, the fastp tool [75] was used to perform adapter trimming and quality control of the sequencing reads. For each quality-controlled set of reads, SPAdes version 3.12 [76] was used for genomes assembly, with C. difficile str. 630 as a conservative reference and the--untrusted-contigs and--careful options. Each assembly was automatically annotated using the Prokka pipeline [77] with a minimum contig length of 200. In order to verify the identity of the assembled genomes as strains of C. difficile , the predicted genes from Prokka were taxonomically annotated using Centrifuge [78] against their pre-compiled index of bacterial, archaeal, viral, and human genomes. Only samples that were clearly identifiable as strains of C. difficile were kept.
To identify the tcdA and tcdB genes from all strains, the phmmer tool was used to search for matches to TcdA (uniprot accession # Q189K5_CLOD6) and TcdB (uinprot accession # Q189K3_CLOD6) as queries. In order to distinguish true sequence variants from poorly assembled, low-quality, or chimeric variants, only hits that clearly represented well-assembled toxin sequences (that is, yielding a protein equal to or greater than 1,800 amino acids in length) were retained. Sequences with apparent N- or C-terminal truncations representing less than 1% of the total assembled data set were also removed. In total, the final re-assembled set of redundant TcdA and TcdB sequences consisted of 3,542 and 5,994 sequences, respectively. Redundancy in each of these sets was removed by clustering with CD-HIT version 4.6 [79] at 100% identity. Non-redundant sets were aligned using the L-INS-i algorithm of the MAFFT package [80].
GenBank homologs of TcdA and TcdB were also identified via a BLAST search of the NCBI non-redundant database on Feb 8, 2020. TcdB and TcdA sequences from C. difficile strain 630 were used as queries. Homologs were filtered to those with E-value < 1e-10, 70% identity and query alignment coverage, which removed partial sequences. In addition, we manually curated 70 reference C. difficile strains collected from previous studies [12,46]. For these 70 genomes, we manually identified corresponding strains within the NCBI or SRA database and identified tcdA and tcdB genes based on pre-computed genome annotations or through similarity searches. Fourteen genomes could not be associated with tcdA and tcdB genes in the NCBI; for these cases, raw genomic reads were retrieved from European Nucleotide Archive (ENA) and were assembled using SPAdes as described earlier.
A combined dataset of TcdA and TcdB homologs was created by pooling SRA-derived, NCBI-nr derived, and the manually curated set of sequences. The combined set of sequences were aligned using MUSCLE [81] with default parameters as implemented in Seaview [82]. Due to significant length variation at the C-terminus of TcdA alignment, only the CROP-less core region (1–1874) of the alignment was kept for subsequent analysis; while the entire TcdB alignment (1–2366) was used. Redundant sequences (100% identity) were removed as well as sequences annotated as partial that contained truncations in the alignment.
TcdA and TcdB alignments were then processed separately using an analysis pipeline implemented within R. For each case, the multiple sequence alignment was converted to a distance matrix using the dist.alignment() function from the seqinr package [83]. Average linkage hierarchical clustering was performed using the hclust() function. Pairwise sequence similarities were mapped onto the clustering tree and visualized using the ComplexHeatmap package [84] and clustering threshold were chosen to generate subtypes with strong internal (within-cluster) and lower external (between-cluster) similarities based on visual analysis and quantitative analysis of percentage identity distributions.
For analysis of amino acid variation, we converted the alignments into data matrices using the alignment2matrix() function from the BALCONY R package [85]. We then identified all variant residues across all alignment positions relative to the sequences of TcdA and TcdB from strain 630 as a reference. Residues implicated in frizzled binding [42], bezlotoxumab binding [63], PA41 binding [66], and E3 binding [65] were then analyzed in terms of their variation across subtypes. E3 binding residues were identified by analysis of PDB structure 6OQ5 [65], by selecting atoms in chain A (TcdB) within a 4 Å distance of chain E (E3) using PyMol’s distance algebra functions. The ComplexHeatmap R package was used for data visualization.
Analysis of repeats in TcdA was done using InterproScan as part of the InterPro 80.0 database [86]. The number of detected matches to ProSite’s cell wall-binding repeat profile (PS51170) was counted in the A1.1 reference sequence (UniProt P16154, 2710 aa), the longest (3070 aa) and shortest (1889 aa) variants of TcdA in our database.
To map sequence conservation on to the structure of TcdB, we used the ConSurf server [87] with the TcdB alignment as input and the recently determined crystal structure of full-length TcdB (PDB ID 6OQ5) [65] as the template. Default parameters (neighbor-joining with ML distance and Bayesian calculation of conservation scores) were used. Structural visualization was done using PyMol version 2.3.4, using the recommended script (https://consurf.tau.ac.il/pyMOL/consurf_new.py) with insufficient data hidden from the image.
We retrieved 2,118 assemblies for 1,934 representative C. difficile genomes (https://www.ncbi.nlm.nih.gov/genome/tree/535) from the NCBI. The snippy pipeline (https://github.com/tseemann/snippy) was used to map all genomes to the reference (strain 630, GCA_000003215). For phylogenetic reconstruction, we analyzed 14,194 SNPs across 88 conserved marker genes (those present in C. difficile ) derived from the PhyEco Firmicutes dataset [88]. A phylogeny was reconstructed using RAxML with the GTR+GAMMA model [89]. All TcdA and TcdB homologs from the NCBI were then subtyped by BLAST against our database of labeled toxin subtype sequences, using only the conserved portion (region 1–1874) of the TcdA alignment, and the full 1–2366 regions from the TcdB alignments. An assignment script written in Perl was used to parse BLAST output files and assign subtypes. The subtype “X” associated with the best matching reference sequence (highest sequence identity) was assigned if the alignment coverage exceeded 90% and labeled as complete; otherwise, it was labeled as a partial sequence.
TcdA and TcdB alignments were analyzed by SplitsTree version 4.0 [57]. A NeighborNet tree visualization was produced using protein maximum-likelihood distances according to the WAG model of evolution. The Phi test for recombination was performed as implemented in SplitsTree which selected a window size of 100 for TcdA with k = 3 and a window size of 100 for TcdB with k = 21.
For visualization of recombinant blocks and haplotype structure within TcdA and TcdB protein alignments, we developed a modified algorithm based on a previous method from Wang et al. [56] for comparative genomic visualization. An implementation of this method in the R programming language is available at https://github.com/doxeylab/haploColor. The algorithms works as follows:
The algorithm was applied directly to the TcdA and TcdB alignments and run for both 4 and 16 iterations (TcdB) and 16 iterations (TcdA).
The DiffBase web server was developed as an R shiny() application. Contained within DiffBase is an implementation of BLAST+. Individual sequences can be submitted to the server, where the blastp program is run to find matches from within the entirety of the server sequence repositories. An E -value cutoff of 1e-10 is used to filter hits, and the results are sorted by percent identity between query and target sequences. Toxin groups can also be viewed in a phylogenetic tree visualized using ggtree R package [90]. Metadata about group members was obtained from the NCBI Identical Protein Group (IPG) database. The source code is freely available at https://github.com/doxeylab/diffBase.
Tweets
Phylogenomics of 8,839 Clostridioides difficile genomes reveals recombination-driven evolution and diversification of toxin A and B
Twitter
Facebook
WhatsApp
LinkedIn
