Dynamic Subspecies Population Structure of Vibrio cholerae in Dhaka, Bangladesh

Cholera has been endemic to the Ganges Delta for centuries. Although the causative agent, Vibrio cholerae, is autochthonous to coastal and brackish water, cholera occurs continually in Dhaka, the inland capital city of Bangladesh which is surrounded by fresh water. Despite the persistence of this problem, little is known about the environmental abundance and distribution of lineages of V. cholerae, the most important being the pandemic generating (PG) lineage consisting mostly of serogroup O1 strains. To understand spatial and temporal dynamics of PG lineage and other lineages belonging to the V. cholerae species in surface water in and around Dhaka City, we used qPCR and high-throughput amplicon sequencing. Seven different freshwater sites across Dhaka were investigated for six consecutive months, and physiochemical parameters were measured in situ. Total abundance of V. cholerae was found to be relatively stable throughout the 6-month sampling period, with 2 × 105 to 4 × 105 genome copies/L at six sites and around 5 × 105 genome copies/L at the site located in the most densely populated part of Dhaka City. PG O1 V. cholerae was present in high abundance during the entire sampling period and composed between 24 and 92% of the total V. cholerae population, only showing occasional but sudden reductions in abundance. In instances where PG O1 lost its dominance, other lineages underwent a rapid expansion while the size of the total V. cholerae population remained almost unchanged. Intraspecies richness of V. cholerae was positively correlated with salinity, conductivity, and total dissolved solids (TDS), while it was negatively correlated with dissolved oxygen (DO) concentration in water. Interestingly, negative correlation was observed specifically between PG O1 and salinity, even though the changes in this variable were minor (0–0.8 ppt). Observations in this study suggest that at the subspecies level, population composition of naturally occurring V. cholerae can be influenced by fluctuations in environmental factors, which can lead to altered competition dynamics among the lineages.


Introduction
Vibrio cholerae is a normal inhabitant of aquatic environments such as rivers, estuaries, and coastal waters and has been detected in diverse geographic locations worldwide [1]. Toxigenic strains of V. cholerae are capable of causing cholera, an acute life-threatening diarrheal disease [2,3], which is a major public health concern because of its high morbidity and mortality with an estimation of 1.3 to 4 million cases and 21,000 to 143,000 deaths worldwide each year [4]. Seven cholera pandemics have struck human civilization so far; the first pandemic of cholera started in 1817 and was followed by six others in the next 200 years, leaving a devastating human death toll [5]. In the environment, V. cholerae

Importance
Although Vibrio cholerae has been studied for decades as a human pathogen, little is known about its abundance and population structure in the environment before an outbreak. In this study, we have employed high-throughput amplicon sequencing of a species-specific region within a protein-coding gene (viuB) to track subspecies lineages of V. cholerae, including the PG lineage responsible for cholera, in an urban region endemic for cholera. Coupled with real-time qPCR, this method provided subspecieslevel resolution of the abundance and lineage composition of V. cholerae populations. This is a key to understanding how intraspecies diversity could play a role in regulating the abundance of pandemic V. cholerae, with environmental factors contributing to the relative fitness of different lineages. is a diverse species, with more than 200 serogroups being identified based on their surface polysaccharide O antigen [6]. However, only O1 serogroup dominates the pandemic generating (PG) lineage [7,8], which has been responsible for all seven cholera pandemics, and other serogroups mostly represent V. cholerae environmental strains which are generally non-toxigenic [5,9]. Serogroup O1 strains of the PG lineage are further classified into two biotypes: the classical biotype that was shown to cause the fifth and sixth pandemics and believed to be associated with the earlier pandemics, and the El Tor biotype which is the causative agent for the seventh pandemic of cholera [5,8].
Cholera remains as an emerging and reemerging disease today [5]. The seventh pandemic of cholera started in 1961 and is still ongoing, known as the world's longest and most persistent pandemic [4,8]. Despite being a global disease, the world's worst cholera epidemics can be traced back to Ganges Delta; hence, that area has been thought to serve as a reservoir for the disease [10,11]. High population density and improper sewage disposal along with uncontrolled industrialization constantly pollute the water sources in this region. Dhaka, the capital city of Bangladesh, is considered as a hyper-endemic region for cholera. The disease continually exists at a low incidence but exhibits biannual seasonal outbreaks in this area [12]. Consequently, Dhaka has been a center for clinical research on cholera for many years. Recent advancement on cholera surveillance intended to increase prevention, preparedness, intervention, and awareness of the disease has been chiefly established based on analysis of clinical cases and clinical samples [13]. However, little attention has been paid to dynamics of the natural population of V. cholerae, which could play an important role in the epidemiological outcome of cholera. As V. cholerae is an autochthonous member of the aquatic environments, understanding its population structure and ecology requires knowledge of the temporal and spatial variations in abundance of various genotypes in its natural habitat.
Most studies looking at the environmental aspects of the disease have been dependent on culture-based isolation of toxigenic V. cholerae [14]. Culture-based studies underestimate V. cholerae abundance and diversity, because of important limitations including enrichment biases and the viable but non-culturable (VBNC) state [15][16][17]. Real-time qPCR analyses of environmental DNA targeting V. cholerae or O1 serogroup-specific genes are quantitative but do not provide information on the subspecies composition of the populations analyzed [18]. Additionally, culture-independent studies, including 16S ribosomal RNA sequencing, can be helpful for genus identification but do not usually provide resolution at the species or subspecies level [19]. A fluorescent antibody staining method has been used to enumerate viable V. cholerae cells in water samples [20], but only targets the O1 or O139 serogroups and cannot distinguish between O1 strains belonging to the PG lineage and strains from other lineages bearing that antigen.
To overcome these limitations, we have applied highthroughput sequencing of a species-specific, highly variable region of a gene, encoding the vibriobactin utilization protein, viuB [21]. This protein releases iron captured by the siderophore vibriobactin inside the cell [22], which is a housekeeping function for V. cholerae. In combination with using qPCR [23], amplification and sequencing of the partial viuB gene allowed for the elucidation of the spatio-temporal abundance of the PG lineage and other lineages belonging to the V. cholerae species in Dhaka's water reservoirs during six consecutive months from October 2015 to March 2016, roughly encompassing after-monsoon, fall, winter, and spring seasons in Bangladesh. The study reveals the continual presence of the PG lineage and occasional reduction of its abundance in conjunction with increases in other lineages, providing insights on the influence of environmental factors in subspecies-level population dynamics of V. cholerae.

Study Site
Surface water samples were collected from seven different locations (site 1 to site 7) in Dhaka (23.8103° N, 90.4125° E), Bangladesh (Fig. 1), biweekly for six consecutive months from October 2015 to March 2016. Dhaka is the capital city of Bangladesh surrounded by a river system mostly composed of four rivers: Turag, Buriganga, Shitalakshya, and Balu. A population of > 21 million was recorded for this area in 2020, with a density of 23,234 people/km 2 within a total area of 300 km 2 ( Fig. 1) (World Population Review, 2020). The physical distance between site 1 and site 7 is shorter (9.9 km) than the distance between site 1 and site 7 to the other five sites (approximately 21 to 25 km). The climate of Dhaka is categorized as tropical wet and dry with a distinct monsoon season. The average water temperature recorded is 26.1 °C (19.1 °C in January and 29.1 °C in June). The visual inspection of the study areas indicated that local markets surrounded these sites, and human intrusion such as bathing, swimming, washing household utensils dishes, as well as bathing domestic animals was frequent. Occasional direct defecation in the water at several study sites was also noticeable.

Sample Collection and Processing
Water samples (200 mL) were collected directly from the sources in sterile Nalgene bottles. Using 50-mL sterile polypropylene syringes, 50 mL water sample was filtered through 0.22-μm Sterivex filters (Millipore). Total DNA extraction from the biomass on the filters was done through the following three consecutive steps: cell lysis and digestion, DNA extraction, and DNA concentrating and washing according to the protocol developed by Wright et al. [24]. Briefly, the filters with trapped biomass are incubated with lysozyme (Thermo Fisher) to lyse the cells and RNase A (New England Biolabs) to remove RNA at 37 °C for 1 h. Afterwards, the filters are treated with proteinase K (New England Biolabs) and 20% SDS at 55 °C for 2 h. Then, the lysate were treated with an equal volume of phenol:chloroform:isoamyl alcohol (IAA) to separate the DNA layer from proteins and debris followed by an equal volume of chloroform:IAA, to remove any remaining phenol from the DNA sample. Final DNA was eluted in TE buffer for downstream applications. To reduce impurities that can act as PCR inhibitors during amplification, all extracted DNA samples were further treated with OneStep PCR Inhibitory Removal Kit (Zymo Research) by following the user manual instructions with 90-180 μL of yield achieved from 100 to 200 μL of the eluted extracted DNA sample. Treated samples were kept at − 20 °C for further analysis.

Culture-Based Identification
Bacterial isolates were isolated as described elsewhere [16]. Briefly, water samples were enriched in alkaline peptone water (APW) (Difco) at 37 °C for 6 to 8 h before plating. About 5 μL of enriched APW broth was streaked, using an inoculating loop, onto thiosulfate-citrate bile salts-sucrose (TCBS) plates at 37 °C for 18 to 24 h. Yellow colonies on TCBS agar were subcultured onto gelatin agar (Difco) plates. Gelatinase-positive colonies were subjected to PCR for amplification of the V. cholerae species-specific ompW gene and O1 serogroup-specific rfbO1 gene for the confirmation.

Physicochemical Parameters
Surface water quality was measured in situ at the sampling sites. EXO2 multiparameter Sonde (YSI, Xylem Brand, USA) allowed for simultaneous measurement of pH, dissolved oxygen (DO), conductivity, total dissolved solids (TDS), salinity, and water temperature. Properly calibrated sensors attached to the instrument were placed in each site while sampling and data were recorded for analysis.  Location of sampling sites in the study. Dhaka, the capital city of Bangladesh, is surrounded by the river system of four rivers including the Turag River, Buriganga River, Shitalakshya River, and Balu River, as indicated by blue pins on the map. Seven different sampling sites are indicated with "S" (from S1 to S7) along with the approximate human population density corresponding to each site (marked by pins with distinct colors). Information on the human population density in this figure was adapted from Khatun et al. [ [25] were used to prepare amplified viuB products for sequencing. Amplified viuB products (2 μL) were used as a template for the tagging PCR reaction with the same reagents as above with a set of forward and reverse primers that contained appropriate Illumina adapters, a sample-specific 8-nucleotide index sequence, a 10-nucleotide pad, a 2-nucleotide linker, and the gene-specific sequence described above, for a total of 70 bp and 65 bp [21]. This tagging PCR reaction was performed as follows: initial denaturation at 98 °C for 30 s, followed by 2 cycles of denaturation at 98 °C for 10 s, annealing at 55 °C for 6 s, and extension at 72 °C for 1 s, and final extension at 72 °C for 1 min. Use of gene-specific primers during amplification and subsequent tagging to create dual-indexed PCR products facilitated improved yield of amplicons and prevented biased amplification due to an unexpected interaction of non-primer sequences with the template. Additionally, eight tagging reactions were done for each sample to obtain an adequate concentration of amplicon DNA for further analysis. All eight reactions of the same sample were pooled together and run on a 2% agarose gel in 1 × Tris-acetate-EDTA buffer where two bands of very similar size, a smaller band (around 360 bp) representing only half-tagged PCR products, and a slightly bigger band (428 bp) of the fully tagged product were visualized. The 428-bp bands were cut out of the gel. PCR products were purified using Wizard SV Gel and PCR Clean-Up System (Promega) according to the instructions by the manufacturer. The concentration of cleaned PCR products was then measured using a Qubit Fluorometer (Thermo Fisher) with a Qubit dsDNA HS Assay Kit (Thermo Fisher) and pooled together in equal concentrations (> 10 ng/μL). The pooled samples were then concentrated using a Wizard SV Gel and PCR Clean-Up System (Promega) according to the instructions by the manufacturer. Quality control of the pooled and the concentrated sample was performed using an Agilent 2100 Bioanalyzer system. Sequencing was performed using a v3 (600 cycles) Illumina sequencing reagent kit and in Illumina MiSeq platform using the sequencing facility at the Molecular Biology Service Unit (MBSU) at University of Alberta, Canada.

Sequence Analysis
Processing of amplicon sequence reads was performed following the procedure described by Kirchberger et al. [21]. Briefly, de-multiplexed raw reads were processed in R [26] using the DADA2 pipeline [27]. Forward and reverse reads were trimmed due to a drop-off in read quality in the first 10 bp as well as after 240 bp and 160 bp for forward and reverse reads, respectively. Assembled overlapping forward and reverse reads were therefore 272 bp in length and 11 bp shorter than the fully sequenced region. Reads with a maximum expected error rate > 1% were also discarded based on DADA2 analysis. After this procedure, 1072 unique sequences remained in the dataset. Chimera detection implemented in DADA2 was then performed on pooled samples, leaving a total of 460 unique sequences. Sequence reads were assigned to viuB alleles, 25 of which were composed of more than 1000 reads (with an average of 100,000 reads per sample) and with the 12 most abundant alleles representing > 99% of all reads were considered for further analysis.

Real-Time qPCR Amplification
A real-time qPCR assay was performed to determine the abundance of total V. cholerae and toxigenic V. cholerae O1 [23]. For total V. cholerae, the viuB gene encoding vibriobactin utilization protein B was targeted using the probe (5′-/56-FAM/TCA TTT GGC/ZEN/CAG AGC ATA AAC CGG T/3IABkFQ/-3′) and forward and reverse primers (5′-TCG GTA TTG TCT AAC GGT AT-3′ and 5′-CGA TTC GTG AGG GTG ATA-3′, respectively), for a 77-bp product. For the V. cholerae O1 serogroup, rfbO1 gene was targeted; probe (5′-/5HEX/AGA AGT GTG/ZEN/TGG GCC AGG TAA AGT/3IABkFQ/-3′), forward primer (5′-GTA AAG CAG GAT GGA AAC ATA TTC -3′) and reverse primer (5′-TGG GCT TAC AAA CTC AAG TAAG-3′) were used for a 113-bp product [23]. For the qPCR reaction, Dynamite qPCR Mastermix which is a proprietary mix, developed and distributed by the MBSU at University of Alberta, Canada, was used. The volume of the PCR reaction was 10 μL containing 5 μL of 2 × Dynamite qPCR Mastermix, 1 μL of each of 500 nM primer-250 nM probe mix, 1 μL of molecular grade water, and 2 μL of DNA template. Realtime quantitative PCR was performed using the following conditions: initial primer activation at 95 °C for 2 min followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min in the Illumina Eco Real-Time PCR System. The qPCR method has been reported to have an analytical detection limit for the viuB and rfbO1 genes as three gene copies per reaction with 1.5 × 10 6 copies/L of water as the lowest detectable number without filtration [23].

Statistical and Multivariate Analyses
Two-dimensional visualizations of the overall V. cholerae community structure were performed using non-metric multidimensional scaling (NMDS) with Bray-Curtis distance at operational taxonomic unit (OTU) levels based on unique viuB alleles. "bioenv" function in the vegan package was used to find a non-parametric monotonic relationship between the dissimilarities in the samples matrix and for plotting the location of each site in a two-dimensional space [28]. Abundance and environmental variable data were analyzed using R. Multiple subplots were generated for the environmental variables against viuB allele richness, and a correlation test was performed to determine any significant relationships between them. Hierarchical clustering analysis was done with the R function hclust to identify clusters of viuB alleles with the most similar distribution based on the abundance data.

Phylogenetic Analysis and Determination of Genome Similarity
For phylogenetic analysis, two representative strains (Table S1) were chosen for each major viuB allele. V. cholerae strains isolated from Dhaka, genome of which has already been sequenced [29], and reference genome sequences were obtained from NCBI database. Wholegenome analysis was done using mugsy v1.2.3 with default parameters [30]. Mugsy outputs were analyzed using the Galaxy web server (Galaxy: a web-based genome analysis tool for experimentalists) (https:// usega laxy. org), and the phylogenetic tree was constructed using RAxML v8.2.11 under the GTRGAMMA model with 100 bootstrap replicates [31]. Defaults were chosen for all other parameters. The phylogenetic tree was visualized using iTOL [32]. In silico digital DNA-DNA hybridization (dDDH), which is a proxy for traditional DNA-DNA hybridization (DDH) value, was used to calculate genomic similarities [33]. All pairwise comparisons for dDDH were calculated using GGDC with default parameters [34]. Allelic differences were also used to determine genomic similarities. To determine allelic differences, each genome was first annotated using RAST [35]. A set of 2443 genes common in most V. cholerae strains were then identified using USE-ARCH [36]. Allele designations and identification were subsequently performed using automated scripts made available by BIGSdb [37]. Finally, pairwise allelic differences for all isolates were calculated using an in-house script and only loci present in both isolates of a pair are considered. Sequences for all gene alleles are available on https:// pubml st. org/ vchol erae under the cgMLST scheme.

The Pandemic Generating Lineage Is the Most Abundant Vibrio cholerae Genotype in Dhaka's Water System
Dhaka, one of the most densely populated cities in the world (> 21 million residents as of June 2020), is located within the Ganges Delta. It is an inland city, with a water system primarily consisting of freshwater rivers and canals. Cholera is endemic there and shows biannual peaks in reported cases during the spring and fall, i.e., before and after the monsoons [38]. However, little is known about the abundance and distribution of V. cholerae in natural waterbodies and if it correlates with environmental and human factors. To track the abundance of total V. cholerae and PG lineage of V. cholerae, we used culture-based detection as well as qPCR analysis. A 272-bp hypervariable stretch of the viuB marker gene, which is present in a single copy in V. cholerae [21], was amplified and sequenced from fortnightly samples taken from seven different water reservoirs in and around Dhaka City from October 2015 to March 2016 (Fig. 1). The presence of V. cholerae could be detected by both culture-based isolation and qPCR throughout the sampling period. Abundance of total V. cholerae did not show large fluctuations with fortnightly sampling, ranging from 2 × 10 5 to 5 × 10 5 genome copies/L (Fig. 2). V. cholerae O1 could only be detected in ~ 8% of the samples using culturebased isolation followed by identification with conventional PCR (Table S2). However, qPCR for the rfbO1 gene using DNA extracted from water detected toxigenic PG V. cholerae O1 throughout the sampling period, with monthly average abundance of this genotype of 1.2 × 10 5 to 2.5 × 10 5 genome copies/L (Fig. 2).
Using viuB amplicon sequencing, 25 viuB alleles were found in total, differing from each other by two or more single nucleotide polymorphisms. Each unique allele represents a specific V. cholerae lineage, roughly the equivalent of a multilocus sequence typing (MLST) clonal complex [21,39]. Among them, 12 viuB alleles (each with > 20,000 sequence reads) were included for further analysis of the V. cholerae community composition, while the other 13 alleles representing < 1% of the total population were excluded following the strategy adopted previously [21]. Total amplicon sequencing reads were normalized based on qPCR data. The count of viuB-73 allele from amplicon sequencing corresponding to the PG lineage [21] significantly correlated with the rfBO1 copy number obtained by qPCR (Pearson correlation P value < 0.05). Therefore, viuB-73 largely corresponds to the 7 th pandemic El Tor V. cholerae in the Dhaka environment and was present in all seven sites throughout the 6 months of sampling (Fig. 2). Whether this allele dominates V. cholerae communities in the water bodies of Dhaka is not surprising, given the endemicity of cholera in the city and the contamination of water bodies with human waste. Interestingly, viuB-73 displayed high abundance based on viuB sequencing throughout the 6-month sampling period, as opposed to only occasional detection of strains corresponding to this genotype using culture-based techniques [15,38]. The sustained presence of PG V. cholerae detected by both qPCR and amplicon sequencing underscores the importance of culture-independent methods of tracking toxigenic V. cholerae in nature. Although regular cholera episodes occur year-round, the bacterium is usually only culturable during the two seasonal peaks [38], presumably being present in its VBNC state at other times [17]. In subspecies-level population analysis, viuB-73 remained the predominant allele in Dhaka, representing 24 to 92% of the total V. cholerae population. In a study of the natural V. cholerae population in cholera-free Oyster Pond on the US east coast, viuB-73 was found sporadically at lower relative abundance and other viuB alleles dominated the V. cholerae population [21]. Even though the current study could not present the picture of a whole year due to our sampling limitations, it revealed that PG V. cholerae was present even in the so-called "non-epidemic" months of the year (December-February) and was the predominant genotype throughout the sampling period in the surface water of cholera-endemic Dhaka. During these months, a majority of the V. cholerae population is thought to remain in VBNC form [15,40], which could be the cause of their disappearance from culture-based surveillance. Another possibility is that there is a lower level of virulence gene expression during this period, due to the lower water temperature. Other environmental or human factors may also lead to the reduction of the number of infections during this period, even though toxigenic V. cholerae survives in the environment. Drastic change in the environment like flooding and excessive rainfall can contaminate the water supply and lead to the rise of clinical cases to pose a threat to the healthcare system.

Subspecies-Level Diversity Correlates with Variation in Environmental Parameters
Studies conducted to look into populations of V. cholerae in natural habitats have shown that different environmental factors might influence their abundance [41,42]. But, how variations in these parameters influence the subspecies-level diversity of V. cholerae is not known. In this study, several environmental variables were measured in situ: pH, DO, conductivity, TDS, salinity, and water temperature. This was to gain insight on their effect on the absolute abundance of the species and the relative abundance of the different viuB alleles found in the water reservoirs in Dhaka. Salinity, pH, and temperature showed little variation at the seven sites sampled, ranging from 0 to 0.8 ppt, from 6.5 to 8.0, and from 30 to 32 °C, respectively (Table S3). Conductivity (134.5-1608 μs/cm) and TDS (67.2-804 mg/L) changed noticeably with sites and the time of sampling (Table S3). DO also changed across the time and space ranging from 0.1 to 5.05 mg/L (Table S3). Potential links between environmental parameters and the diversity (richness) and population composition were assessed using correlation analyses. Diversity (richness) of V. cholerae alleles was significantly correlated with salinity, conductivity, TDS, and DO (Fig. 3,  Fig. S1). Salinity, conductivity, and TDS were found to be positively correlated with allele richness (Pearson correlation coefficients 0.44, 0.41, and 0.40, respectively; P < 0.01), while DO (mg/L) showed a negative correlation (Pearson correlation − 0.30; P < 0.05) with that population characteristic. In previous studies, various physicochemical variables have been found to be associated with the abundance and persistence of V. cholerae in aquatic environments and with the risk of cholera outbreaks [43,44]. Temperature and salinity have been observed to influence planktonic populations, which is a well-known habitat for V. cholerae in the aquatic ecosystem [41]. Adaptation to a wide range of salinity levels also facilitates V. cholerae's survival in various aquatic environments (from coastal to inland water) [42,43]. Previous studies suggest that abundance of V. cholerae decreases with increasing salinity and that they are most abundant in salinities ranging from 0 to 10 ppt [42,43]. Multiple linear mixed effect regression analysis was performed to find out which of the independent parameters impacted the richness of the subspecies population significantly. A repeated measure analysis was performed using a site-level random intercept term. When all the six environmental variables were used in the model, it revealed that salinity, conductivity, and TDS were highly correlated with each other. In the final model, only four parameters (pH, DO, salinity, and WT) were therefore considered as independent predictors of the variance in viuB allele richness ( Table 1). The only variable that was found to be significant in that regards (at significance level of 0.05) is salinity (P value = < 0.001) ( Table 1), increasing richness by 19.1% for each increase of 0.1 ppt. As conductivity and TDS were strongly correlated with salinity (Pearson's rho = 0.93 for both), the effect of those two parameters on richness can be expected to be similar to salinity.
Redundancy analysis (RDA) was done using the Monte Carlo permutation test with the aim of testing the significance of the constraint ranking model to reveal what environmental factor/s impacted the shift in viuB allele composition. In the RDA biplot, two axes explained 27.1% of the variation in allele composition (Fig. 4). Salinity was the strongest factor that was positively correlated with allele richness. DO was negatively correlated, and temperature only showed moderate positive correlations with viuB allele richness. viuB-05, viuB-06, and viuB-07 were strongly correlated with each other. The distribution and abundance of the cluster containing these three alleles was positively correlated with salinity and negatively correlated with viuB-73. This latter allele had a clear negative correlation with salinity as observed in the RDA plot, whereas the distribution and abundance of viuB-05, viuB-06, and viuB-07 were positively correlated with salinity.
In a study of the V. cholerae population in cholera-free Oyster Pond coastal ecosystem (Falmouth, MA, USA), the viuB-73 allele was rarely found in the ocean [21], whereas it was present in the brackish pond and lagoon water connected to the ocean, suggesting that high salinity might represent an environmental barrier to the dispersal and range of cholera. The negative correlation between viuB-73 abundance and salinity observed in this study also suggest that PG V. cholerae (viuB-73) might be adapted to low salinity. The PG lineage represented by viuB-73 might have higher tolerance of rapid variations at low salinities, as it was consistently predominating in all but one site of Dhaka where salinity levels fluctuated between 0 and 0.8 ppt. The only site at which viuB-73 was not dominant was site 7 (Fig. 5), in which salinity was more stable and only varied from 0.4 to 0.5 ppt and never dropped below 0.4 ppt. This suggests that even at the subspecies level, small environmental fluctuations could have an influence on the population composition of V. cholerae. In the low-lying Ganges Delta, salinity intrusion is considered a major threat due to the changing climate. Reduced upstream discharge, sea level rise, and other catastrophic events such as cyclones can lead to an increase in the salinity of inland water bodies [43]. A salinity increase of about 26% was recorded in coastal regions in Bangladesh over the last 35 years [45]. Such a shift in salinity could be affecting the composition of V. cholerae populations, possibly changing the distribution and abundance of various lineages, some of which could pose a threat to human health.
DO, a measure of free non-compound oxygen present in an aquatic system, is another variable influencing microbial communities [46]. A study in Hood Canal (Washington, USA) demonstrated that there was a strong negative correlation between bacterial richness and DO [47]. So far, there are no studies describing the correlation between the intraspecies diversity of V. cholerae and DO, but it was demonstrated earlier that V. cholerae is most abundant in low-DO environments [48,49]. In our analysis, DO was found to impact the V. cholerae diversity negatively in the Dhaka environment (Fig. 3, Fig. 4). DO could also be one of the factors contributing to the differential adaptation at the subspecies level. Among the sampling locations, site 7 had the lowest average DO (0.90 mg/L). Average DO in site 1 was also low at 0.95 mg/L, but the other five sites (site 2-6) had noticeably higher average DO concentration ranging from 2.7 to 3.0 mg/L. Conductivity, salinity, and TDS usually have a negative relationship with DO, which might have caused the lower DO observed in sites 1 and 7. Site 7 was the only one where viuB-73 was outcompeted in abundance by other viuB alleles (viuB-05 and viuB-06) (Fig. 5). V. cholerae, being a facultative anaerobe, can adapt to the low-oxygen conditions by utilizing alternative energy-producing pathways (i.e., nitrate utilization) [50]. It is possible that some lineages are at an advantage in lower oxygen concentrations, giving them the ability to co-exist with other lineages that usually outcompete them, leading to increased diversity at lower DO.  Fig. 4 Redundancy analysis (RDA) illustrating the relationships between Vibrio cholerae community at the subspecies level and environmental variables. Environmental factors and viuB allele richness were marked by red arrows relative to the first two canonical axes (RDA1 and RDA2). The angles in the biplot between viuB alleles and environmental variables, and between viuB alleles themselves or envi-ronmental variables themselves, reflect their correlations. An angle less than 90° suggests positive correlations, and an angle approaching 180° suggests a strong negative correlation between variables. The length of the arrow measures the degree of effect an environmental variable has on the community, and a longer arrow indicates a greater effect on the community Temperature is another factor affected by climate change that is known to have an impact on V. cholerae populations. While this species is found at a wide range of temperatures (10 to 30 °C), the highest abundance is observed at > 20 °C [42]. However, no significant correlation between temperature and viuB allele diversity (Pearson correlation coefficient 0.1) was observed in Dhaka (Fig. 3). Water temperatures remained within the range of 27.4 to 30.8 °C throughout the sampling period, a pattern that differs dramatically from regions where shifts in temperature are noticeable during summer and winter [51]. In a temperate region, attachment to particles and hosts has been shown to increase when temperature increases above 22 °C, contributing to changes in the lineage composition of a V. cholerae population when a seasonal change occurs in a temperate climate [21]. It is likely that at conditions with consistent high temperatures such as those found in a tropical region like Bangladesh, V. cholerae is not overly responsive to this parameter, either in terms of its growth rate or particle attachment behavior. Overall, this study unravels the role of important environmental parameters on the composition of the V. cholerae community in Dhaka. It must be noted that observations from this study are based on biweekly sampling at a single point per site. We have also lost some samples due to difficulties in transportation, which led to a lack of continuity in temporal data (Fig. 5). As a result, this study is likely to have missed small fluctuations in environmental determinants and population composition. A high-resolution time series analysis would pinpoint on how even small daily fluctuations in those parameters impact the V. cholerae community at a subspecies level and the possible ecological outcome of those changes in the population composition.  viuB gene copies/L

High Human Population Density Correlates with Changes in the V. cholerae Population of Dhaka Reservoirs
Dhaka is one of the most densely populated areas in the world, with a density of 23,234 people per km 2 within a total area of 300 km 2 . This huge burden of human population has significant impact on the ecology and evolution of V. cholerae and, consequently, on the epidemiology of cholera in this region [52,53]. Based on the demographic records of Dhaka City, the area surrounding sampling site 7 (Kamrangirchar) is the most densely populated (100,000 people/km 2 ) among the seven sites studied, whereas population density at the other six sites ranged from 10,000 to 60,000/km 2 (Fig. 1) [54]. Thus, human impacts on the water reservoir in this area is expected to be much higher than other sites, with a higher level of fecal contamination and industrial waste mostly from tannery industrial units. This area is mostly inhabited by a dense low-income population where people use shared hygiene facilities such as showers and toilets [55]. Open defecation has been reported from poorly maintained shared facilities. Additionally, high population density in this area frequently causes an overload of septic tanks, which results in the overflow of untreated effluent to the water reservoir [56]. The open drainage system commonly causes mixing between sewage and fresh water, increasing the possibility of V. cholerae transmission between the water reservoir and local population. Among the seven different sites studied, most of them had a TDS concentration of < 300 mg/L on average, except sites 1 and 7, where TDS varied from 308 to 472 mg/L and from 476 to 575 mg/L, respectively (Table S3). TDS value comes from the combination of the disassociated electrolytes and other compounds such as dissolved organic matter [57]. Typically, natural bodies of water have dissolved solids due to the dissolution and weathering of rocks and soil. However, human activities also influence the concentration of TDS in water reservoirs, which is especially likely in Dhaka's inland water bodies because of urban runoff as well as wastewater discharge. This heavy influence of human population at Kamrangirchar likely led to different population dynamics of V. cholerae at this site compared to other locations in the city. However, in our sampling sites, conductivity, salinity, and TDS were found to be tightly correlated, indicating inorganic ions might be the major contributor of the variation in TDS values. Hence, major cause of higher conductivity, salinity, and TDS values observed in site 7 can also be indicators of increased chemical pollution there in comparison to other sites. Noticeable differences in the diversity and population composition in site 7 might suggest that human population density can have a direct or indirect impact on the V. cholerae population at the subspecies level. Although absolute abundance of total V. cholerae remained stable at all seven sites throughout the 6-month sampling period, there are some for site 7, where total V. cholerae abundance was around 28 to 43% higher than that in any other site (Fig. 2). It is unclear if this higher abundance is due to a more constant input of V. cholerae from human waste, an indirect increase in numbers from nutrients linked to this waste, or the existence of a niche allowing the expansion of a particular lineage that does not compete with others (sympatry). Unlike sites 1 to 6, at which viuB-73 was more abundant than all other alleles combined (Fig. 5), three viuB alleles (viuB-05, viuB-06, and viuB-07) had greater combined abundance than viuB-73 at site 7 (Fig. 5). The trend is most pronounced during December 2015 to March 2016, coinciding with a decrease of water quality in the reservoir during December to April [58]. These lineages have so far not been found in environmental surveys outside of Dhaka. Interestingly, these three alleles (viuB-05, viuB-06, and viuB-07) have been found in strains representing a basal long-branch clade in the whole-genome phylogeny of known global V. cholerae strains [29]. Representative isolates of this clade were found to be indistinguishable from typical V. cholerae based on conventional phenotypic tests, but they are phylogenetically and genotypically divergent [59].
Allelic differences and in silico DDH of whole-genome sequences show that organisms represented by these three alleles (viuB-05, viuB-06, and viuB-07) are not more closely related to each other on average than most pairs of V. cholerae strains would be (Fig. S2). Despite this lack of genetic similarity, their presence and absence are strongly correlated, as they are usually co-occurring (Fig. 6, Fig. 4). Surprising dominance of the viuB-05, viuB-06, and viuB-07 lineages over viuB-73 exclusively at the Kamrangirchar (site 7) location suggests that this group of lineages is specifically adapted to the environmental conditions found at this site. As the site also stood out from a human population density point of view, and that these lineages were most abundant in this site, it suggests a potential human link to the ecology of these lineages in Dhaka. Indeed, strains phylogenetically related with these lineages have been isolated from human samples in different parts of the world [59]. Abundance of those lineages was positively correlated with TDS, conductivity, and salinity, which can be considered as indicators of the impact of human intrusion in the water. Among the seven sampling sites, these three related parameters were highest on average at site 7. These three alleles (viuB-05, viuB-06, viuB-07) were most abundant at site 7 (54% of total V. cholerae), where viuB-73 had the lowest abundance among all sites (24%). This suggests that some form of competition could be taking place between different V. cholerae lineages and/or that the lineages respond differently to environmental factors present at this site. Recently, in a more robust study of genomic, phylogenetic, and phenotypic characterization of isolates harboring these alleles has suggested for it to be a novel sister species of V. cholerae and named as Vibrio paracholerae. Possible link of the human population to the abundance and distribution of this novel species makes this lineage of bacteria compelling candidates for future studies looking into the ecology of human-adapted genotypes of V. cholerae.

Intraspecies Interaction Could Influence Relative Abundance of PG V. cholerae O1
Even though overall viuB-73 was the predominant allele over the 6-month sampling period, spatial and temporal analyses indicate that direct or indirect intraspecies competition among V. cholerae genotypes could play an important role in the population dynamics observed in Dhaka. Whenever viuB-73 displayed a drop in abundance, another lineage carrying a viuB-05, viuB-06, viuB-07, viuB-25, viuB-45, viuB-51, or viuB-79 allele underwent a rapid expansion. Although the abundance of these other alleles was not directly quantified by qPCR, the absolute abundance measurement of the total V. cholerae and PG V. cholerae O1 by qPCR (Fig. 2) supported the observation coming from amplicon sequencing data. Furthermore, this shift in abundances is repeated multiple times within our dataset. For example, viuB-25 increased in relative abundance when the viuB-73 abundance decreased at sites 1 and 7 (Fig. 5). Notably, a new viuB allele, viuB-79 (for which no cultured isolates have been found), also appeared when viuB-73 abundance was reduced at sites 6 and 7 (Fig. 5). All of these alleles showed an inverse correlation with viuB-73 relative abundance (Fig. 6). Hierarchical clustering analysis of viuB amplicon sequencing data (Fig. 6) shows that statistically significant clustering (Hopkins statistics (H) < 0.5) only occurred between viuB alleles representing strains harboring viuB-05, viuB-06, and viuB-07 alleles (H < 0.5). viuB-73 showed a significant negative correlation with most other viuB alleles (H > 0.5) [60]. The only viuB alleles positively correlated with viuB-73 were viuB-10 and viuB-39 (Fig. 6), the latter being a ubiquitous allele present at low abundance across sampling sites as well as other geographical locations (Fig. 5).
The cause of these shifts in the abundance of some alleles is unclear. It could be due to a differential response to environmental factors or trophic interactions. In the Oyster Pond Fig. 6 Co-occurrence of V. cholerae lineages in Dhaka water reservoirs. The Hopkins statistics (H) was used to assess the clustering tendency of various viuB allele sequences in the dataset. The threshold value was 0.5, meaning that if H < 0.5, data showed significant clustering. Visualization of analyzed viuB amplicon sequence data is presented with blue and brown color gradients (H < 0.5, indicating that the data is highly clusterable which is shown with a brown color gradient, and H > 0.5, indicating that the data is not clusterable which is described by blue color gradient). Significant clustering (H < 0.5) is observed between viuB alleles 05, 06, and 07. Single asterisk indicates significant clustering (MA, USA) V. cholerae population, divergent responses of different lineages were observed in ocean, lagoon, and pond, where ecosystem parameters varied substantially [21]. It is plausible that V. cholerae lineages responded differently to even the slight environmental variations observed at the sites sampled in Dhaka. Differential response of lineages to phage predation can also be a cause of shifts in the community composition. Effect of predation by bacteriophages on the V. cholerae population composition can be modulated by the environmental factors, i.e., nutrient availability [61], which, in turn, can give advantage to certain lineages to outcompete others under certain environmental conditions. Another factor influencing this differential response could be the ability of various strains to avoid predation. The abundance of V. cholerae is influenced by grazing from heterotrophic protists [41]. To overcome the grazing pressure, V. cholerae executes different strategies such as morphological shift, i.e., from smooth to rugose, resulting in the production of Vibrio polysaccharide (VPS) that helps to encase themselves in biofilm and resist predation [61]. V. cholerae can also survive predation by becoming intracellular in a range of amoeba [62]. They are also able to kill grazers using T6SS [63]. This system encodes a syringe-like structure that can pierce cellular envelopes of other bacteria and some eukaryotes, injecting effector proteins that can kill the recipient if it does not process the cognate immunity protein [64]. This phenomenon could influence the population composition of V. cholerae, with different lineages having varying effectors and immunity proteins providing varying predatory success [64,65]. Another possibility is that the incompatibility between subspecies likely plays a role in the diversity and dynamic of the V. cholerae populations in Dhaka. Most of these genotypes (T6SS effector and immunity protein profiles) suggest they are incompatible and can kill each other on contact [66]. Even though the PG lineage carrying viuB-73 could seemingly outcompete strains represented by all other alleles at sites 1-6, it was outcompeted by the three co-occurring alleles corresponding to the viuB-05, viuB-06, and viuB-07 group at site 7. This site differs from the other six sites in terms of environmental parameters and surrounding human population density. These observations suggest that environmental conditions can play an important role in shaping the intraspecies competition of V. cholerae in their natural environment to impact diversity-and subspecies-level population dynamics of the V. cholerae.

Conclusion
For centuries, the Ganges Delta has been a reservoir for pandemic and non-pandemic V. cholerae lineages. Cholera endemicity is usual in this area, and Dhaka is one of the most densely populated megacities in the world, with cholera epidemics occurring biannually before and after the rainy season. This study revealed the consistent presence of PG V. cholerae lineage in the natural water bodies of Dhaka at a considerable proportion of the total V. cholerae population, suggesting that the toxigenic bacterium is circulating year-round in this urban aquatic environment or is constantly shed by human carriers. Moreover, population analysis with subspecies-level resolution revealed that other V. cholerae lineages co-existed with PG V. cholerae O1 in that environment. Intraspecies niche specialization and potentially subspecies interactions with these other lineages could decrease PG V. cholerae O1 abundance occasionally, especially when environmental parameters, such as consistently higher TDS, salinity, or lower dissolved oxygen concentration, favored other lineages. Variability of these parameters, even on a small scale, correlated with changes in V. cholerae population composition and diversity, with more temporal stability of salinity and TDS at one particular site correlating with a dramatically different lineage composition. The same site also displayed the highest human population density among the seven sampling locations, indicating humans as a possible contributing source of organic and inorganic materials, or of V. cholerae themselves, altering V. cholerae population composition. Consistent human interaction and frequent leaking of sewage into the water being linked to a substantially different lineage composition would suggest that the human gut may also serve as a potential reservoir for PG O1 and other lineages of V. cholerae, resulting in year-long persistence of V. cholerae belonging to these lineages through the transmission cycle between human and aquatic environments. A study of the microbiomes of individuals living near water reservoirs is essential to identify if humans play a direct or indirect role in the differences observed between the aquatic V. cholerae populations found at different Dhaka sites.
Technology-Futures and Faculty of Graduate Studies and Research to TN, MTI, and KYHL.
Data Availability All the genome sequences and relevant genomic and epidemiological data related to the isolates used are publicly available on PubMLST (https:// pubml st. org/ vchol erae/).