Genomic attributes of *Vibrio cholerae* O1 responsible for 2022 massive cholera outbreak in Bangladesh

Vibrio cholerae, native to aquatic environments and the causative agent of cholera, has evolved continually across regions. The current seventh pandemic El Tor (7PET) lineage originated from non-pathogenic ancestors, acquiring tcpA, CTX prophage, VSP-1, and VSP-II, enabling global spread since 1961. The Ganges Delta is a historical hotspot for pandemic clone evolution. Genomic studies show 7PET spread from the Bay of Bengal in at least three overlapping waves, with numerous transcontinental transmissions. Multiple transmission lineages have been defined in Africa and Latin America, and two contemporary 7PET wave-3 lineages circulate in Asia. In Bangladesh, two lineages, BD-1 and BD-2, have been identified with significant genomic differences. In 2022, Dhaka experienced its largest cholera outbreak in 20–25 years, motivating genomic investigation of the causative strains to understand evolutionary dynamics and drivers of the outbreak.

Prior work defined multiple waves of global 7PET transmission and numerous lineages across continents, including 13 transmission lineages in Africa and three in Latin America. In South Asia, Bangladesh and India have reported co-circulating 7PET wave-3 lineages. Comparative genomics highlighted differences among lineages in SNP/indel burden, mobile genetic elements (e.g., SXT ICE conferring multidrug resistance), VSP islands implicated in environmental responses, and phage defense elements like PLE. Earlier Bangladeshi studies (1991–2017) identified BD-1 and BD-2, with BD-2 showing greater SNP/indel accumulation, enriched AMR genes, mobile cassettes, and PLEs, and predominance in endemic cholera during 2013–2017.

Clinical sampling and sequencing: 21 V. cholerae O1 El Tor isolates were collected from patients at icddr,b, Dhaka, between February and April 2022. These were compared with 267 genomes from the authors’ collection (1991–2021) and 693 publicly available genomes (1957–2017), totaling 981 genomes for certain analyses. Whole-genome sequencing: Genomic DNA was extracted using Qiagen kits; QC and quantification used NanoDrop 1000 and Qubit 4.0. Libraries (300–350 ng DNA) were prepared with Illumina DNA Prep and sequenced as 150 bp paired-end reads on an Illumina NextSeq 500. Bioinformatics: Raw reads were quality filtered with fastp. Assemblies were generated using SPAdes v1.5.4.3 and scaffolds with Ragout v2.3. Annotation was performed with Prokka v1.14.5. AMR gene profiles were identified using ResFinder and ABRicate v0.8.1. Core non-repetitive, non-recombinant SNPs (n=6399) across 981 genomes informed maximum likelihood phylogenies; an outgroup (V. cholerae N16961) was used for rooting. Temporal phylogenomics employed BEAST v2.6.7 (MCMC). Lineage/group-associated SNPs were detected via genome-wide association using Pearson’s chi-squared tests with contingency tables; analyses were implemented with in-house R scripts. rfb gene sequences were extracted by BLAST (using references JX565645–JX565867). Functional enrichment used PANZER and GSEA-Pro v3.0; GO term networks were built with REViGO. Phenotypic testing: Additional 2022 isolates (March–September) were tested for cexB allele via DMAMA-PCR and for antimicrobial susceptibility via disk diffusion assays. Phylogenetic and lineage definitions: Maximum-likelihood and unrooted trees delineated Bangladeshi subclades BD-1, BD-1.1, BD-1.2 within a global clade, and BD-2 within an Asian clade. Comparative analyses included African transmission lineages (e.g., T12–T13) and Latin American lineage LAT-3.

- The 2022 Dhaka outbreak strains clustered within the globally distributed 7PET wave-3 clade and formed a distinct Bangladeshi subclade, BD-1.2, separate from earlier BD-1 and BD-1.1 and from the previously dominant BD-2 Asian clade. - Temporal dynamics: Negative correlation in predominance between BD-1 and BD-2 from 2001–2007; co-detection 2007–2015; BD-2 predominated 2013–2017; BD-1.2 surged from 2018 onward, becoming responsible for the 2022 outbreak. - Phylogenomics: Analysis of 981 genomes and 6399 high-quality core SNPs showed BD-1.2 closely related to Indian strains (circa 2016) and to six defined lineages; BEAST dated the MRCA of BD-1 and BD-1.2 around 2005 (95% HPD). - Genomic distinctions: BD-1.2 and BD-2 differ in SXT ICE types, VSP-I variants, VPI-1 status, gyrA/parC alleles, presence of PLE, and ctxB carriage. BD-1.2 outbreak strains harbored ctxB, while BD-2 often carried PLE (phage defense) absent in most BD-1.2. - Clade associations: BD-2 clustered within an expanded Asian clade; BD-1.2 clustered within a global clade with strains from multiple regions (including T11–T13, LAT-3). - AMR profiles: Global clade lineages (including BD-1.2) showed a largely uniform resistance gene pattern distinct from Asian clade lineages. In the Asian clade, tetA(D) occurred in 71–99% of strains; dfrA13 and gyrA(C) in 17–46%; floR in 20% of AS-2 but not in descendant IND-2/BD-2. All lineages carried parC Ser82 mutations; BD-1 and BD-1.2 had different resistance alleles. A BD-2 strain exemplified extreme multidrug resistance to 11 tested drugs. - Serotype switching: BD-1.2 isolates were Ogawa (2016–2019) but switched to Inaba predominance in Bangladesh from September 2020 through 2022. Sequence analysis linked Inaba to insertions at rfb (e.g., 20 bp or 12 bp insertions at positions 27–28). - Epidemiology: BD-2 was associated with endemic cholera (especially 2013–2017), whereas BD-1.2 is newly recognized and associated with the massive 2022 Dhaka outbreak. - Additional testing: All additional 2022 isolates tested carried cexB and exhibited tetracycline resistance consistent with BD-1.2 profiles.

The study addressed why Bangladesh experienced a massive cholera outbreak in 2022 by demonstrating that a globally distributed 7PET wave-3 subclade, BD-1.2, expanded locally and displaced the previously dominant BD-2 lineage. This lineage shift likely altered the epidemiology due to BD-1.2’s distinct genomic features (e.g., ctxB presence, different SXT/VSP/VPI configurations, absence of PLE in most BD-1.2) and a contemporaneous serotype switch from Ogawa to Inaba. The serotype switch may have increased the susceptible population due to reduced cross-protection between serotypes, potentially contributing to higher case counts. Temporal phylogenomic analysis indicates BD-1.2’s MRCA emerged around 2005, with subsequent spread linked to strains in India, consistent with transboundary transmission dynamics observed in prior 7PET waves. AMR analysis further differentiates global and Asian clades, with BD-2 harboring extensive resistance and PLE, while BD-1.2 aligns with a global clade resistance profile. Altogether, these findings link evolutionary shifts, serotype dynamics, and genomic determinants to the resurgence and scale of cholera in Bangladesh in 2022.

This work identifies BD-1.2, a subclade within the global 7PET wave-3 lineage, as the principal driver of the 2022 cholera outbreak in Dhaka, Bangladesh, displacing the previously dominant BD-2 lineage. The study delineates key genomic and epidemiological differences between BD-1.2 and BD-2, including AMR gene content, presence/absence of PLE, SXT/VSP/VPI configurations, and a serotype switch from Ogawa to Inaba linked to rfb insertions. Phylogenomic analyses suggest a recent common ancestry for BD-1/BD-1.2 around 2005 and highlight lineage turnover over time. Future research should monitor ongoing lineage dynamics, investigate mechanisms and impacts of serotype switching on population immunity and transmission, and assess the role of mobile genetic elements and phage interactions in shaping outbreak potential and antimicrobial resistance.