Clostridioides difficile infection (CDI) is an urgent public health threat with limited therapeutic options. In this work, we developed a messenger RNA (mRNA)-lipid nanoparticle (LNP) vaccine targeting C. difficile toxins and virulence factors. This multivalent vaccine elicited robust and long-lived systemic and mucosal antigen-specific humoral and cellular immune responses across animal models and independent of changes to the intestinal microbiota. Vaccination protected mice from lethal CDI in both primary and recurrent infection models, and inclusion of non-toxin cellular and spore antigens improved decolonization of toxigenic C. difficile from the gastrointestinal tract. Our studies demonstrate mRNA-LNP vaccine technology as a promising platform for the development of novel C. difficile therapeutics with potential for limiting acute disease and promoting bacterial decolonization.

mRNA design and production

The amino acid sequence of the TcdA and TcdB receptor binding domain (RBD) containing the combined repetitive oligopeptides (CROP) motifs were obtained from GenBank accession numbers P16154.2 and P18177.3 respectively. The amino acid sequences of PPEP-1 and CdeM was obtained from GenBank accession number Q183R7, and WP_009893169.1, respectively.  The putative N-glycosylation sites were disrupted by substituting the asparagine residue at the predicted N-glycolysation site with a glutamine (N to Q) to prevent posttranslational epitopes masking in eukaryotic cells. The sequences underwent codon optimization and GC enrichment using our proprietary algorithm to improve expression and reduce potential immunogenicity of the in vitro transcribed mRNA. The codon-optimized sequences were gene synthesized by Genscript with an optimized and modified IL-2 secretion signal, and cloned into our proprietary in vitro transcription template containing an optimized T7 promoter, 3’UTR, 5’UTR and a 100-adenine tail. The TcdA, TcdB, PPEP-1, and CdeM nucleoside modified mRNA sequences were prepared using the MegaScript transcription kit (ThermoFisher Scientific), co-transcriptionally capped using the CleanCap™ system (TriLink Biotechnologies) and purified using a modified cellulose based chromatography method (36), precipitated, eluted in nuclease free water, and quantified using the NanoDrop One system. Length and integrity were determined using the Agilent BioAnalyzer 2100 system. Endotoxin content was measured using the GenScript Toxisensor chromogenic assay, and values were below detection levels (0.1 EU/mL). mRNA was frozen at -20oC until formulation.

Production and characterization of mRNA-LNP vaccines

The hydrodynamic size, polydispersity index (PDI) and zeta potential of mRNA-LNPs were measured using a Zetasizer Nano ZS90 (Malvern Instruments, Malvern, UK). The mRNA encapsulation efficiency was determined using a modified Quant-iT RiboGreen RNA assay (Invitrogen). Endotoxin levels were determined using the Limulus Amebocyte Lysate (LAL) chromogenic assay and found to be <0.5 endotoxin unit (EU)/mL.

Phylogenomic analyses

Genomic analyses were conducted on 137 C. difficile isolates sourced from the work of Bushman, Frederic D., et al(12). Initial quality control was implemented with FastQC v0.12.1 and MultiQC v1.14 (37). Subsequently, de novo assembly was performed on the quality-assessed raw reads utilizing Shovill v1.1.0 (38). The genome assemblies were subjected to quality evaluation with CheckM v1.1.6 and BUSCO v5.4.7(39, 40). Genome annotations were then performed employing bakta v1.8.1(41). The core-genome alignment was carried out using Roary v3.13.0(42),  which subsequently was used to infer the maximum likelihood phylogeny in RAxML-NG v1.1(43) using the general time-reversible (GTR) substitution model (44) accounting for among-site rate heterogeneity using the Γ distribution and four rate categories(45) for 10 individual searches with maximum parsimony random-addition starting trees and random topology starting tree. Node support was evaluated with 100 nonparametric bootstrap pseudoreplicates (46). The phylogeny was rooted at the midpoint using the Tree Of Life (iTOL) v6.8 interactive tool(47). For subsequent visualization and annotation of the phylogeny, ggplot2, and ggtree were utilized (48, 49).

Comparative Genomic Analysis of Protein Domains

The amino acid sequences for the TcdA CROP domain, TcdB CROP domain, PPEP-1, and CdeM (excluding the IL-2 secretion signal domain) were extracted from UniProt Knowledgebase with primary accession numbers of P16154 from position 1840 to position 2710, P18177 from position 1851 to position 2366, Q183R7 from position 27 to position 220, WP_009893169.1 from position 2 to position 164, respectively. These sequences served as the seed sequences for constructing a protein-type BLAST database. Subsequently, BlastX was applied to execute a comparative analysis of the amino acid sequences in the database against the 137 genomes included in the phylogenetic assessment. The hits with the highest bit score from BlastX were utilized to determine the percentage of coverage and identity(50-57).

Identification of Toxigenic Strains via Comparative Genomic Analysis

The coding sequences for tcdA, tcdB, tcdC, tcdE, tcdR, cdu1, and cdd1 were sourced from the comprehensively annotated complete genome of C. difficile 630 (NCBI Reference Sequence: NC_009089.1). The nucleotide sequence for the PaLoc region was inferred and extracted by the locations of the flank genes cdu1 and cdd1. These nucleotide sequences served to construct an original nucleotide-specific BLAST database. To enhance the likelihood of detecting diverse toxin genes, additional sequences of toxin genes and the PaLoc region were extracted from 54 representative toxigenic strains and incorporated into the database. The raw sequencing reads from 54 representative toxigenic strains were downloaded and assembled (58). The initial round of blastn was conducted against the genomes of these 54 toxigenic strains using the original database, setting the parameters to accommodate more gaps and mismatches to boost the coverage of hits. The sequences from the highest-scoring hits were subsequently extracted and used as seed sequences in a new BLAST database. Genomes of good quality were also used to infer and extract sequences from the PaLoc region using flanking genes cdd1 and cdu1. A conserved non-coding region of 115 bp was employed as a marker for non-toxigenic strains (59). A second round of blastn was performed against the 137 genomes included in the previous phylogenetic analysis, using sequences from the new BLAST database. An arbitrary threshold of 80% coverage and 80% identity was applied to the highest-scoring hits and used as criteria for determining the presence or absence of genes. Strains were classified as non-toxigenic if they exhibited the 115 bp conserved region and lacked toxin genes or the PaLoc region, and vice versa.

In vitro transfection of mRNA-LNP and Western blot

Neuro2a cells were seeded in 12-well plates at a density of 0.2x106 cells per well and transfected 24 hours post-seeding at a confluence of around 80%. mRNA-LNPs (1 mg/mL) were diluted (1:10) in PBS and added to each well at a final dose of 2.5 µg/well. Medium was aspirated, and cells were directly lysed 24 hours post-transfection using the EZLys tissue protein extraction reagent (BioVision #8002-500) supplemented with 1X protease inhibitors (Sigma-Aldrich, #0493116001). Untreated cells were used as controls. Lysate was clarified at 15,000 rpm for 5 minutes using a refrigerated centrifuge to pellet cell debris, and protein concentration in the supernatant was determined using the Pierce microBCA assay. A total of 12 μg of cell lysate was loaded into a 4-15% polyacrylamide gel (Bio-Rad #4561083) and transferred to a PVDF membrane using the ThermoFisher iBlot system (dry transfer). Membranes were blocked with 5% skimmed milk extract, washed, and incubated for 2 hours with primary antibodies against TcdA and TcdB (Abcam 19953 and 252712) at a final concentration of 4 µg/mL. PPEP-1 was detected using a polyclonal serum from immunized mice (Boosted sera at 1:100 dilution, incubation 4 hours). Membranes were then washed three times using 1X TBST and incubated for 1 hour in the presence of an HRP-conjugated donkey anti-goat IgG (Abcam 97040) to detect TcdA and TcdB. PPEP-1 was detected using mouse polyclonal sera. Proteins were visualized using ECL detection reagent (Cytiva #RPN2209) on a GE ImageQuant™ system.

Cell viability measurements

Cells that are relevant for toxicity assessment and representative of different tissues that come in to contact with the expressed immunogens encoded on our mRNA-LNPs were used for in vitro transfection and assessment of viability post-transfection. HUVEC (CRL-1730), Caco-2 (HTB-37), and human primary skeletal muscle cells (PCS-950-010) were cultured as per ATCC using F12K Medium (ATCC 30-2004), EMEM (ATCC 30-2003), and Mesenchymal Stem Cell Basal medium (ATCC PCS-500-030) respectively. The muscle cells were cultured in the presence of the Primary Skeletal Muscle Growth Kit (ATCC PCS-950-040). Cells were seeded in 96-well plates at a density of 10,000 cells per well and transfected with 1 or 3 µg TcdA or TcdB mRNA-LNP.  Full-length recombinant TcdA (8619-GT-020) or TcdB (6246-GT-020), recombinant Luciferase mRNA-LNP (Luc LNP), and untreated cells were used as controls. After incubation for 24 or 48 hours, 50 µL XTT cell proliferation reagent (Roche 11465015001) was added to each well, and plates were incubated for an additional 5 hours before measuring absorbance at 492 nm using the Thermofisher Varioskan plate reader. Cell viability was calculated relative to the untreated cells.

Generation of recombinant TcdA, TcdB, PPEP-1, and CdeM proteins

The region encoding the TcdA and TcdB receptor binding domains (including CROP domain), the full length PPEP-1 coding sequence, and the full length CdeM coding sequence was codon-optimized for expression in E. coli, gene synthesized at Genscript, and cloned into the pET30a vector. An N-terminal 6XHIS tag followed by the TEV recognition/cleavage site was introduced after the start (ATG) codon of all three constructs to allow for purification using affinity chromatography. E. coli BL21(DE3) was transformed, and positive clones were used to inoculate 1 L cultures. Recombinant proteins from the supernatant were concentrated and purified using the NI-IDA column on an AKTA Avant 150 system. Purified proteins were cleaved with TEV protease, buffer-exchanged into PBS and 5% sucrose (pH 7.4), filter-sterilized, quantified using the micro-BCA assay (Pierce), and aliquoted and stored at -80ºC. Aliquots were tested for precipitation following multiple freeze and thaw cycles and the purity was determined using densitometric analysis of a Coomassie blue stained SDS-PAGE gel under reducing conditions. All proteins displayed more than 90% purity and had the predicted molecular weights using Western blot.

Spleen and lymph node harvests

Spleens were collected, processed as single cells, filtered using 70 µm cell strainers in complete RPMI 1640, centrifuged, and red blood cells lysed in ACK lysis buffer to obtain a clear single cell suspension. Splenocytes were resuspended in 1 mL complete RPMI 1640 media, counted, and used immediately for studies. Draining lymph nodes (inguinal and popliteal) were collected and processed as described above (without ACK). Cells from the dLNs were resuspended in 250µL RPMI 1640 media, counted, and used immediately.

Production of fluorescently labeled recombinant proteins for antigen-specific B cells

Fluorescently-labeled recombinant TcdA, TcdB, and PPEP-1 protein were prepared using the Lightning-Link R-Phycoerythrin (R-PE) and Lightning-Link (R) Rapid Alexa Fluor 647 conjugation reagents (Novus Biologicals, 703-0010 and 336-0005). TcdA, TcdB, and PPEP-1 were diluted to 0.5 mg/mL in PBS and 50μg reacted in the presence of 1:10 (v/v) of LL modifier for 3 hours at room temperature. The labeling reaction was stopped in the presence of 1:10 (v/v) LL quencher for 30 minutes and stored at 4°C until use.

 Flow cytometry analysis of T and B cells

Tfh cells: 2 million cells were stained with anti-mouse CD16/32 antibody for 20 minutes and stained with an anti-CXCR5- biotin for 30 minutes on ice, washed twice and incubated with Streptavidin BV-421 in the presence of surface antibodies (Table S5) for 30 minutes, washed, fixed and permeabilized with the FoxP3/Transcription Factor Staining Kit (eBioScience) according to the manufacturer instructions and stained for Bcl6. Following intracellular staining, cells were washed twice, fixed in 300 µL 1% paraformaldehyde for acquisition. The gating strategy, as well as the antibody list and catalog numbers are provided (Table S5; fig. S6A).

Germinal Center (GC) and Memory B cells: 2 million cells per sample incubated with anti-mouse CD16/32 antibody for 20 min at 4°C. Cells were then washed with FACS buffer (2% FBS in PBS) and stained for 1 h using antibodies (Table S3). Following staining, cells were washed twice, and fixed in 300 µL 1% paraformaldehyde for acquisition. The gating strategy, as well as the antibody list, fluorescent TcdA and TcdB RBD probes, and catalog numbers are provided (Table S3, fig. S6B-D).

T cells: 2 million splenocytes were stimulated with 2.5 µg/mL of TcdA CROP, or TcdB CROP, or PPEP-1 peptide pools (15 mers, 4 amino acid overlapping peptide pool) in a FACS tube for 6 hours at 37oC, 5% CO2 with 2 mg/mL anti-CD28 (Tonbo, 40-0281-M001) providing co-stimulation. Stimulations proceeded for 1 hour before adding 5 mg/mL brefeldin A (Biolegend, 420601), 2 mM monensin (Biolegend, 420701), and 5 mg/mL anti-CD107a (Biolegend, 121610) Alexa Fluor 647 for 5 hours. DMSO served as a negative control and the combination of 50 mg/mL phorbol 12-myristate 13-acetate and 1 mg/mL ionomycin served as a positive control. After a total of 6 hours, samples were washed with PBS, stained with Live/Dead Aqua for 5 minutes, blocked using anti-mouse CD16/32 antibody for 20 minutes, and stained extracellularly for 30 minutes using antibodies (Table S4). Cells were washed in FACS buffer, fixed and permeabilized using the Cytofix/Cytoperm kit (BD Biosciences, 554714), and stained intracellularly using antibodies for 30 min (Table S4). Following intracellular staining, cells were washed twice, and fixed in 300 µL 1% paraformaldehyde for acquisition. The gating strategy, as well as the antibody list and catalog numbers are provided (Table S4; fig. S9A).

All samples were acquired on a BD LSR II equipped with 4 laser lines and 18 PMTs and data were analyzed in FlowJo v10.

In vivo toxin neutralization assays

To determine the LD100 of intraperitoneally injected TcdA and TcdB, recombinant toxins were prepared at doses of 5, 10, 25, 50, 100, and 125ng in 200μL of sterile PBS.  Naïve 7-week-old C57BL/6J mice were randomly assigned into groups of 5. Recombinant proteins were injected i.p. and mice were monitored for behavior, body condition, and mortality every 4 hours post-injection. Doses were selected for use in subsequent studies by their ability to kill 100% of animals within 4 hours after injection. To determine neutralization capacity of vaccine-elicited immune response, mice were immunized twice i.m. with 1 or 5mg of TcdA or TcdB monovalent mRNA-LNP or monovalent recombinant protein vaccine with or without alum adjuvant (details on vaccine preparation below). Two weeks after last immunization, immunized mice were either challenged i.p. with five times the previously determined LD100 (625ng rTcdA, 125ng rTcdB) and monitored for behavior, body condition, and mortality every 4 hours post-injection or terminally bled. Serum from immunized mice was diluted 1:20 in PBS and incubated with rTcdA or rTcdB LD100 (150 or 25 ng respectively) for 1 hour at 37 °C. Naïve mice were then challenged i.p. with recombinant toxin alone, monovalent immunized serum + recombinant toxin, or control (unvaccinated) serum + recombinant toxin and monitored for behavior, body condition, and mortality every 4 hours post-injection.

 Mouse immunization and C. difficile infection

5-week-old male and female C57BL/6J mice (Jackson Laboratories, strain no. 000664) were used in most studies and BALB/c (Jackson Laboratories, strain no. 000651) were used where indicated. All mice were maintained under specific pathogen-free conditions at either the Children's Hospital of Philadelphia or the University of Pennsylvania. Multivalent mRNA-LNPs (1-5 mg) were mixed at 1:1 w/w ratio, diluted up to 50μL in PBS, and administered intra-muscularly into the hind leg within two hours of thaw. Recombinant protein vaccines were mixed 1:1 with alhydrogel (aluminum hydroxide 2% w/v). 28 days later, mice were either euthanized for analysis of vaccine-induced humoral and cellular immune responses after single immunization or they were boosted with the same dose. 14 days later, mice were either euthanized for analysis of vaccine-induced immune responses after second immunization or they were infected with C. difficile as previously described (28). Antibiotic treatment was administered by providing 0.5g/L cefoperazone in their drinking water ad libitum for 5 days, followed by a 2-day recovery period before C. difficile infection via oral gavage. Two different C. difficile strains were utilized where indicated: VPI 10463 at a dose of 1x106 spores and CD196 at a dose of 1x105 spores. Mice were monitored daily for survival and weight loss, and mice were euthanized when weight loss exceeded 20% of their original body weight.

For studies to analyze immune responses vaccination after C. difficile infection, 5-week-old C57BL/6J mice were infected with 1x105 spores C. difficile CD196 as described above. 14 days post-infection, mice were immunized with 1mg mRNA-LNPs. 28 days later, mice received a second immunization. 14 days later, serum and feces were collected for analysis of vaccine-induced humoral immune responses.

For long-term studies, C57BL/6J mice were immunized as described above. 40 days after last immunization, mice were either euthanized for analysis of vaccine-induced immune responses or infected with C. difficile as described above.

For re-challenge studies, C57BL/6J mice were re-infected with C. difficile 200+ days after primary infection as described above.

Hamster studies

Syrian hamsters (HsdHan: AURA, Envigo, catalog no. 8901M) were maintained at University of Texas Medical Branch. 4-5-week-old male golden Syrian hamsters were pre-bled and vaccinated intra-muscularly with trivalent mRNA-LNP vaccine (1:1:1 w/w ratio), trivalent recombinant protein + alhydrogel, or empty lipid nanoparticles at week 0 and 3. Body weight and clinical signs were collected at days 0, 1, 3, 5, 7, and 14 post-prime and at day 1, 2, 3, 7, and 14 post-boost. Feces were collected 14 days after each immunization. Serum was collected at days 0, 21, 28, and 45 post-prime. Hamsters were sacrificed as per university protocols and observed for gross abnormalities – no findings were recorded.

Detection of antibodies in serum and feces

96-well High Bind StripwellTM plates (Corning) were coated overnight with 1 mg/mL of purified TcdA (CROP), TcdB (CROP), PPEP-1, or CdeM. Plates were washed once with wash buffer (0.5% Tween-20 in PBS) and blocked for two hours at room temperature using a solution of heat-inactivated, IgG-depleted, protease-free bovine serum albumin (2% w/v BSA in PBS). After blocking, plates were washed three times with wash buffer. Feces were homogenized in 1 mL PBS and centrifuged for 10 minutes at 10,000 x g. Sera or fecal supernatants were serially diluted in the blocking solution, added to plates, and incubated for two hours at room temperature. Plates were washed three times before the addition of horseradish peroxidase-conjugated anti-mouse secondary antibody specific to total mouse IgG (1:10,000) and IgA (1:5,000), anti-rhesus IgG (1:20,000), or anti-hamster IgG (1:8,000) in blocking buffer. Plates were incubated for 1.5 hours and washed three times before the addition of tetramethylbenzidine (TMB) substrate solution. The reaction was stopped by adding 2 N sulfuric acid, and the absorbance was measured at 450 nm using a SpectraMaxTM 190 microplate reader. Antigen-specific antibody end-point dilution titer was defined as the highest dilution of serum or feces to give an OD greater than the cut-off OD value determined using the Frey Method (60). Fecal antibody titers were normalized to gram of feces.

C. difficile enumeration and toxin titers

C. difficile burdens were quantified by collecting fecal samples at indicated timepoints and plating on taurocholate cycloserine cefoxitin fructose agar (TCCFA). C. difficile toxin titers were quantified in the stool using a previously described Vero cell cytotoxicity assay (28). Briefly, fecal samples were homogenized in sterile PBS, pelleted, and supernatant was filtered through a 0.2µm filter. Supernatant was diluted serially and incubated overnight with Vero monolayers (ATCC-CCL-81). Toxin titers in stool were calculated as the reciprocal value of the highest dilution that rounded 100% of the cells and normalized per gram of feces.

Histology and scoring

Samples were fixed in 10% neutral buffered formalin (NBF) prefilled HistoTainer™ II (Simport, M96160FW), dehydrated in graded ethanol series, cleared with Xylene, and embedded in paraffin. Sections (5 µm) were collected on Superfrost™ Plus stain slides (Fisher Scientific, Ottawa, ON, Canada), and stained with Hematoxylin and Eosin. Slides were scanned using a NanoZoomer digital slide scanner (Hamamatsu, Boston, MA, USA) and visualized using the NDP® view 2.0 software (Hamamatsu, Boston, MA, USA). Sections were scored in a blind manner by a pathologist based on previously described criteria (30). Histological scores were reported as a cumulative score of three independent criteria: inflammation, edema, and epithelial cell damage.

DNA extraction from stool

Mice were co-housed for one week before experimental manipulation. Stool samples were collected at indicated timepoints for microbiota analysis. Microbial genomic DNA was extracted using DNeasy Power Soil Kit (Qiagen) according to manufacturer’s instructions.

16S rRNA gene library prep 

Barcoded PCR primers annealing to the V4 region of the 16S rRNA gene were used for library generation.  PCR reactions were carried out in duplicate using Q5 High-Fidelity DNA Polymerase (NEB, Ipswich, MA).  Each PCR reaction contained 0.5μM of each primer, 0.34 U Q5 Pol, 1X Buffer, 0.2 mM dNTPs, and 5.0μl DNA in a total volume of 50μl.  Cycling conditions were as follows: 1 cycle of 98ºC for 1 minute; 20 cycles of 98ºC for 10 seconds, 56ºC for 20 seconds, and 72ºC for 20 seconds; and 1 cycle of 72ºC for 8 minutes.  After amplification, duplicate PCR reactions were pooled and then purified using a 1:1 volume of SPRI beads.  DNA in each sample was then quantified using PicoGreen and pooled in equal molar amounts.  The resulting library was sequenced on the Illumina MiSeq using 2x250 bp chemistry.  Extraction blanks and DNA-free water were subjected to the same amplification and purification procedure to allow for empirical assessment of environmental and reagent contamination. Positive controls, consisting of five artificial 16S gene fragments synthesized in gene blocks and combined in known abundances, were also included.

Bioinformatics processing and statistical analysis

Sequence data were processed using QIIME2 (61). Read pairs were processed to identify amplicon sequence variants with DADA2 (62).  Taxonomic assignments were generated by comparison to the Silva reference database version 132 (63), using the naïve Bayes classifier implemented in scikit-bio (64). A phylogenetic tree was inferred from the sequence data using MAFFT (65).  Similarity between samples were assessed by weighted and unweighted UniFrac distance (66, 67). Data files from QIIME were analyzed in R environment for statistical computing. Linear mixed effects models were used at each time point to estimate the mean difference between study groups. Cage information was added to the models as the random effect to account for the coprophagic nature of mice. Relative abundances of bacteria were log10 transformed and modeled as the outcome. Only the bacteria with at least 1% mean relative abundance across samples were tested. Community-level differences between sample groups were assessed using the PERMANOVA test (24, 68). When multiple tests were done, p-values were corrected for false discovery rate using Benjamini-Hochberg method. 

Non-human primate studies

An 18-year-old male rhesus macaque (N = 1) was injected i.m. in the right deltoid with 200 μg of tetravalent mRNA-LNP vaccine (TcdA, TcdB, PEPP-1, and CdeM) in 250μL total volume (1:1:1:1 w/w or 50µg/immunogen). A second age-matched male rhesus monkey (N = 1) was i.m. injected in the right deltoid with 200 μg of an eDHFR-tagged SARS-CoV-2 S2P mRNA-LNP vaccine. NHPs received two immunizations separated by 21 days, and blood was collected at day 1, 8, 21, and 35 post injection. Plasma and PBMCs were separated and stored until use.

Statistical analysis

Data are presented as mean ± SD or mean ± SEM where indicated. Graph schematics and illustrations were created using iTOL v6.8 (47), SPICE 6 (69) and Graphpad prism V10.0.2

References

36.       M. Baiersdorfer et al., A Facile Method for the Removal of dsRNA Contaminant from In Vitro-Transcribed mRNA. Mol Ther Nucleic Acids 15, 26-35 (2019).

37.       P. Ewels, M. Magnusson, S. Lundin, M. Kaller, MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047-3048 (2016).

38.       Seemann, T. Shovill v 1.1.0. Zenodo. https://doi.org/10.5281/zenodo.13351686. (2016).

39.       M. Manni, M. R. Berkeley, M. Seppey, F. A. Simao, E. M. Zdobnov, BUSCO Update: Novel and Streamlined Workflows along with Broader and Deeper Phylogenetic Coverage for Scoring of Eukaryotic, Prokaryotic, and Viral Genomes. Mol Biol Evol 38, 4647-4654 (2021).

40.       D. H. Parks, M. Imelfort, C. T. Skennerton, P. Hugenholtz, G. W. Tyson, CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res 25, 1043-1055 (2015).

41.       O. Schwengers et al., Bakta: rapid and standardized annotation of bacterial genomes via alignment-free sequence identification. Microb Genom 7, 000685 (2021).

42.       A. J. Page et al., Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics 31, 3691-3693 (2015).

43.       A. M. Kozlov, D. Darriba, T. Flouri, B. Morel, A. Stamatakis, RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics 35, 4453-4455 (2019).

44.       C. Lanave, G. Preparata, C. Saccone, G. Serio, A new method for calculating evolutionary substitution rates. J Mol Evol 20, 86-93 (1984).

45.       Z. Yang, Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J Mol Evol 39, 306-314 (1994).

46.       J. Felsenstein, Confidence Limits on Phylogenies: An Approach Using the Bootstrap. Evolution 39, 783-791 (1985).

47.       I. Letunic, P. Bork, Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res 49, W293-W296 (2021).

48.       G. Yu, Using ggtree to Visualize Data on Tree-Like Structures. Curr Protoc Bioinformatics 69, e96 (2020).

49.       G. Yu, T. T. Lam, H. Zhu, Y. Guan, Two Methods for Mapping and Visualizing Associated Data on Phylogeny Using Ggtree. Mol Biol Evol 35, 3041-3043 (2018).

50.       S. F. Altschul, W. Gish, W. Miller, E. W. Myers, D. J. Lipman, Basic local alignment search tool. J Mol Biol 215, 403-410 (1990).

51.       S. F. Altschul et al., Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25, 3389-3402 (1997).

52.       G. M. Boratyn et al., Domain enhanced lookup time accelerated BLAST. Biol Direct 7, 12 (2012).

53.       G. M. Boratyn, J. Thierry-Mieg, D. Thierry-Mieg, B. Busby, T. L. Madden, Magic-BLAST, an accurate RNA-seq aligner for long and short reads. BMC Bioinformatics 20, 405 (2019).

54.       C. Camacho, G. M. Boratyn, V. Joukov, R. Vera Alvarez, T. L. Madden, ElasticBLAST: accelerating sequence search via cloud computing. BMC Bioinformatics 24, 117 (2023).

55.       C. Camacho et al., BLAST+: architecture and applications. BMC Bioinformatics 10, 421 (2009).

56.       A. Morgulis et al., Database indexing for production MegaBLAST searches. Bioinformatics 24, 1757-1764 (2008).

57.       Z. Zhang, S. Schwartz, L. Wagner, W. Miller, A greedy algorithm for aligning DNA sequences. J Comput Biol 7, 203-214 (2000).

58.       K. E. Dingle et al., Evolutionary history of the Clostridium difficile pathogenicity locus. Genome Biol Evol 6, 36-52 (2014).

59.       V. Braun, T. Hundsberger, P. Leukel, M. Sauerborn, C. von Eichel-Streiber, Definition of the single integration site of the pathogenicity locus in Clostridium difficile. Gene 181, 29-38 (1996).

60.       A. Frey, J. Di Canzio, D. Zurakowski, A statistically defined endpoint titer determination method for immunoassays. J Immunol Methods 221, 35-41 (1998).

61.       E. Bolyen et al., Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol 37, 852-857 (2019).

62.       B. J. Callahan et al., DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13, 581-583 (2016).

63.       C. Quast et al., The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41, D590-596 (2013).

64.       N. A. Bokulich et al., Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2's q2-feature-classifier plugin. Microbiome 6, 90 (2018).

65.       K. Katoh, D. M. Standley, MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol 30, 772-780 (2013).

66.       C. Lozupone, R. Knight, UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol 71, 8228-8235 (2005).

67.       C. A. Lozupone, M. Hamady, S. T. Kelley, R. Knight, Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities. Appl Environ Microbiol 73, 1576-1585 (2007).

68.       M. J. Anderson, A new method for non-parametric multivariate analysis of variance. Austral Ecology 26, 32-46. (2001).

69.       M. Roederer, J. L. Nozzi, M. C. Nason, SPICE: exploration and analysis of post-cytometric complex multivariate datasets. Cytometry A 79, 167-174 (2011).