Aeromonas in South Asia: genomic insights into an environmental pathogen and reservoir of antimicrobial resistance

Domain: neuroscience
Published: 2026-01-31
Editorial outcome: minor_revision

Abstract
Aeromonads are an ecologically versatile group of bacteria that cause infections in aquatic animals and are recognised as emerging human pathogens. Despite this, our understanding of Aeromonas diversity, especially the relationship between clinical and environmental strains, remains limited. Here, we present a genomic analysis of the Aeromonas genus, comprising 1853 genomes, and a detailed comparison of clinical and environmental strains from South Asia, including 996 newly sequenced genomes from Bangladesh and India. Phylogenetic analyses revealed that Aeromonas is a highly diverse genus, with no distinct clade separating clinical and environmental isolates. We identified 28 Aeromonas species and 905 novel sequence types, comprising 72.5% of the genomes. Notably, we show a high incidence of antimicrobial resistance (AMR) genes across all isolates, including against front and last-line antibiotics. Finally, we highlight frequent misidentification of Aeromonas as Vibrio cholerae, which is relevant to cholera-endemic regions where both genera co-exist and are associated with diarrhoeal disease. Our study underscores Aeromonas as an important environmental AMR reservoir and emerging multi-species pathogen capable of spilling over into human populations.

Introduction
Aeromonads (within the Aeromonadaceae family) are Gram-negative, facultatively anaerobic bacteria with oxidase and catalase activity. Aeromonas spp. are abundant in aquatic habitats, with over 30 species identified1. At least 19 of these species are linked to human infections and are capable of causing significant diarrhoeal outbreaks. Aeromonas infections are common in children under 5 years old, especially in low- and middle-income countries (LMICs), where there is inadequate access to water and sanitation, resulting in high environmental exposure. The most common symptoms of Aeromonas infection include diarrhoea (which can be mistaken for cholera2,3), localised soft tissue infection, and bacteraemia4–8. Bacteraemia is most common in patients who have underlying health conditions, such as hepatobiliary disease, cancer, or diabetes9. There appear to be geographic differences in prevalence, with <1% of cases of moderate to severe diarrhoea (MSD), attributable to Aeromonas in four African sites and Kolkata, India, whereas Aeromonas was implicated in >22% of MSD cases in children from Karachi, Pakistan, and Mirzapur in Bangladesh10. Aeromonas spp. comprise important pathogens of fish and aquatic animals, causing diseases such as haemorrhagic septicaemia, furunculosis, and motile Aeromonas septicaemia (MAS) in species like salmonids, catfish, and tilapia, often leading to high mortality and significant economic consequences for aquaculture11–17. Additionally, Aeromonas spp. have also been isolated from terrestrial animals and are recognised as zoonotic, with human infections commonly occurring through contact with infected animals or contaminated water, especially via wounds18–20.

Pathogenicity in Aeromonads is complex; they possess a variety of virulence factors that contribute to biofilm formation, cell adherence, invasion, and cytotoxicity. These include polar and lateral flagella21,22, adhesins23, lipopolysaccharides, iron-binding systems24,25, and numerous extracellular toxins and enzymes26 secreted by various systems such as type 2 and type 3 secretion systems27,28. Additionally, quorum sensing systems play a crucial role in colonisation and disease development29–31.

Notwithstanding the above, Aeromonas spp. are frequently misidentified as Vibrio cholerae in clinical and environmental samples. This is explained by sharing many biochemical and microbiological properties with other genera such as Aerobacter, Pseudomonas, Escherichia, Proteus, and Vibrio3,32. Both Vibrio and Aeromonas genera are oxidase-positive, able to thrive in aquatic and marine systems and exhibit similar colony morphologies on Thiosulphate Citrate Bile Salt Sucrose (TCBS) medium routinely used to select for V. cholerae33. These misidentifications poses serious epidemiological challenges, especially in regions where V. cholerae is endemic. Misidentification cannot only lead to misreporting of cholera cases, but also a failure to implement necessary public health control measures. Hence, it is important to consider Aeromonas spp. in cholera surveillance and control efforts in endemic regions.

Accurate identification of Aeromonas spp. remains challenging due to overlapping phenotypic traits and the limited specificity of conventional biochemical tests and selective media, which often produce inconclusive results34–36. Molecular methods such as 16S rRNA sequencing are helpful but complicated by intragenomic variability, reducing species-level accuracy37,38. Emerging technologies like MALDI-TOF MS offer rapid and more reliable genus- and species-level identification37,39,40, while whole-genome sequencing provides definitive species-level resolution and accurate subtyping. Currently, a tiered diagnostic approach combining initial MALDI-TOF screening with phenotypic assays and, when necessary, molecular confirmation at reference laboratories is recommended to ensure accurate detection and effective epidemiological surveillance of Aeromonas spp.40. It is likely that, as for many other pathogens, genomic approaches will replace these techniques where accessible and in settings with well-developed supply chains.

Aeromonads are known to exhibit high levels of antimicrobial resistance (AMR), especially against β-lactam antibiotics. This resistance is primarily attributed to chromosomally encoded inducible β-lactamases in species such as A. hydrophila, A. sobria, and A. salmonicida, conferring resistance to antibiotics including ampicillin, bacitracin, cefoxitin, and piperacillin/tazobactam41–45. Recent genome data from 447 Aeromonas isolates from Pakistan revealed that human-associated Aeromonas in both clinical cases and age-matched controls were highly diverse and displayed a high level of AMR, particularly extended-spectrum beta-lactamases (ESBLs)46. Moreover, Aeromonas serve as natural reservoirs for mobile colistin resistance genes (mcr)47,48, raising concerns about their role in spreading resistance to last-resort antibiotics like colistin across bacterial populations in aquatic and clinical environments. Looking into determinants that may explain the difference in clinical presentations, this study looked across an array of known or putative virulence factors, of which only two, maf2 and lafT, were weakly associated with the MSD cases. However, by focusing on isolates from humans, this study was unable to establish whether human pathogenic strains were a specialised subset of those found in the environment or if Aeromonas spp. infections are truly opportunistic and linked to environmental exposure46.

Given the ubiquitous nature of Aeromonads in the environment, their potentially emerging role in MSD in South Asia, and the fact that they are often cultured from samples taken from patients with suspected cholera, we created a baseline population genomics snapshot of Aeromonas species from clinical and environmental samples. We sequenced two collections of Aeromonas isolates from V. cholerae endemic regions of Bangladesh (2004–2016), and Northern India (2020–2023). These isolates were obtained from environmental sources—including drinking water, ponds, rivers, lakes and drains, as well as from human faeces, together comprising 996 novel Aeromonas genomes. We compared their genetic diversity, predicted virulence gene profiles, and AMR profiles, and contextualised this with genomes from a global Aeromonas population dataset.

Results
To examine the genomic diversity of Aeromonas strains across clinical and environmental settings in South Asia, we included the largest and only available clinical genome datasets from Pakistan and environmental datasets from India and Bangladesh (this study). To place our study into a global context, we also incorporated genomes representing broad geographic and species diversity within the Aeromonas genus. In total, we assembled 1853 Aeromonas genomes, including 996 newly sequenced genomes, of which 198 (10.68%) were mistakenly collected as Vibrio cholerae—133 from suspected cholera cases/households contacts in the ChoBI7 programme in Dhaka city49 and 65 from an environmental surveillance study for V. cholerae conducted in Bangladesh. Further, we generated 798 Aeromonas genomes from isolates taken from 274 water samples collected in Northern India. We contextualised our dataset using public archives, including 441 genomes of clinical origin from children with moderate to severe diarrhoea and corresponding healthy controls, taken in Pakistan46, as well as 416 globally distributed50 Aeromonas genomes. The latter were derived from a wide range of sample types, including human stool, urine, wounds, as well as aquatic environments (e.g., rivers and wastewater), and animals (e.g., fish, duck, dog). The largest proportion of publicly available genomes originated from Denmark (n=100), followed by the United States (n=71) and China (n=62) (see Supplementary Table 1 for detailed information).

We used GTDB-Tk (Genome Taxonomy Database) to perform species assignment, identifying 28 Aeromonas spp. across the 1853 genomes, of which eight species dominated. A. caviae (n=608) was the most numerous in our collection, followed by A. veronii (n=566), A. dhakensis (n=186), A. salmonicida (n=176), A. hydrophila (n=124), A. enteropelogenes (n=105), A. jandaei (n=25), A. sanarellii (n=16), with the remaining 47 Aeromonas genomes representing 20 other Aeromonas species (Fig. 1; summarised in Supplementary Fig. 1a). To understand how our species predictions mapped to established species we used Average Nucleotide Identity (ANI) to evaluate their genetic relatedness. The intra-species ANI values ranged from 94.06% to 100%, falling slightly below the established definition for a single species (>95% ANI), in some instances. Conversely, inter-species ANI values ranged from 79.66% to 95.51% with higher-than-expected ANI and lower levels of expected divergence between A. bestiarum and A. piscicola (95.14–95.51%) (Supplementary Fig. 2).

Looking at the geographic distribution of Aeromonas species across both clinical and environmental samples from India, Pakistan and Bangladesh, where we had the densest sampling (1438 Aeromonas genomes; Supplementary Fig. 1b), Aeromonas spp. found across all three South Asian countries included A. caviae, A. dhakensis, A. enteropelogenes, A. jandaei, A. veronii and A. sanarellii, with A. sanarellii and A. jandaei being rare in all three countries. Whilst overall there were similarities in the relative abundance of the different species, despite differences in sampling depth and strategy (see methods), there were also clear differences. In particular, A. caviae predominated in Bangladesh and Pakistan (31.81%, 63/198 and 64.39%, 284/441; respectively), whereas A. veronii dominated the samples collected in India (48.31%, 386/799 of isolates). Others species, such as A. schubertii, were exclusively seen in Bangladesh (n=1), while A. media (n=4) and A. rivipollensis (n=5) were exclusively seen in India (Supplementary Fig. 1b). What is more, 129 environmental samples from Northern India included between 1–5 different Aeromonas species (Supplementary Fig. 3a). For countries where we had less dense sampling, A. hydrophila was the most common in the United States (42.25%, 30/71), whereas A. salmonicida was the dominant species in the samples from Denmark (99%, 99/100) and China (96.77%, 60/62), (Supplementary Table 1).

In India, all samples were collected as part of a cholera surveillance study. Samples were taken from known cholera hotspots and at times of the year when cholera cases were more common. All water samples were cultured on V. cholerae selective media (TCBS; see methods). Using MALDI-TOF, it was possible to look at the Aeromonas spp. positivity rate, as well as for the presence of a range of other enteric pathogens (Supplementary Fig. 4). Out of 408 collected water samples from Northern India, 274 samples (67.15%) tested positive for Aeromonas spp., 145 (35.53%) for V. cholerae, 112 (27.45%) for E. coli, 198 (48.52%) for Klebsiella spp. and 35 (8.57%) for non-cholera Vibrio spp. (8.57%, 35/408; including species V. fluvialis, V. navarrensis, V. cidicii, V. injensis, and V. metschnikovii).

Of the samples positive for Aeromonas spp., V. cholerae was also isolated from 91 samples (33.57%, 91/274). These co-occurrences were primarily found in river samples (75/91), followed by ponds (13/91) and canals (4/91), across different regions of Northern India. A total of seven Aeromonas species: A. veronii, A. caviae, A. hydrophila, A. dhakensis, A. jandaei, A. enteropelogenes, and A. sanarellii – were identified as cohabiting with V. cholerae. However, there was no clear association between the co-occurrence of V. cholerae with any one specific Aeromonas species; these data largely mirrored the relative isolation rates for Aeromonas species from Northern India (Supplementary Fig. 3b; Supplementary Fig. 3c).

To infer the phylogeny of the genus, we standardised the annotation of our underlying dataset using PROKKA and then defined the Aeromonas spp. pangenome using Panaroo. This comprised 47,542 coding sequences (CDS), of which 1969 were defined as core (found in more than 99% of genomes) and 118 as soft core (found in 95–99% of genomes), with 45,455 accessory genes (found in 0–95% of genomes). We constructed a maximum likelihood phylogenetic tree from a core gene alignment (Fig. 1; see methods). We used Bayesian Analysis of Population Structure (BAPS, implemented in fastBAPS) to delineate genomic clusters and correlated these with ANI values and the species assignments from GTDB-Tk. From BAPS, whilst most taxa fell on deeply branching nodes in the Aeromonas spp. tree with high concordance with the ANI-defined species blocks (Fig. 1; Supplementary Fig. 2), A. veronii fell across two BAPS clusters (15 and 16). Conversely, BAPS clusters 1, 2, 3, and 7 each span multiple species of Aeromonas (Supplementary Table 2). These inconsistencies may reflect differences between the way the BAPS model clusters genomes, which, unlike ANI, is not distance-based.

Next, to determine how much of the diversity within known taxonomic groups had been previously observed, we constructed in silico MLST51 profiles. Of the 1853 MLST profiles generated, only 495 (26.71%) genomes were assigned to 224 known MLST allelic profiles, with 1343 genomes representing 905 novel sequence types (STs). After assignment by pubmlst.org (https://pubmlst.org/organisms/aeromonas-spp), 1838 of our Aeromonas spp. genomes were assigned to 1129 STs (15 genomes remained unassigned; see Supplementary Data 1), now accessible through the Aeromonas MLST database52. Of note, these novel STs defined discrete branches within all species, while 836 STs represented singletons, with clusters of genomes belonging to the same ST being rare. The exceptions to this were ST-2 (n=103) and the newly defined ST-1799 (n=60) for A. salmonicida, as well as ST-2398 (n=20) and ST-251 (n=19), representing A. caviae and A. hydrophila genomes, respectively.

Several virulence genes have been previously described amongst Aeromonas species, including the haemolytic toxins: aerolysin-related cytotoxic enterotoxin (act)41,53,54, heat-labile cytotonic enterotoxin (alt)53–55, heat-stable cytotonic toxins (ast)53,55, hemolysin (hlyA)56, cytotoxin aerolysin (aerA)53,55,56 as well as polar flagellum (fla)55,57,58, lateral flagella (laf)59,60, the zinc-and-iron-dependent metalloendopeptidase elastase (ela)57, the lipase (lip)41,55 and the type III secretion system (TTSS) encoded by ascF-G57. To determine whether any of the known virulence genes are linked to the particular Aeromonas species, we used in silico PCR to understand the distribution of these genes across 1853 Aeromonas species genomes.

We analysed patterns of virulence gene distribution using principal component analysis (PCA) and hierarchical clustering (Fig. 2a; Fig. 2b). The PCA revealed that certain species, such as A. hydrophila, A. dhakensis, A. salmonicida, and A. veronii, formed relatively well-defined clusters, suggesting distinct virulence profiles (Fig. 1; Fig. 2a). Some overlap was observed between the phylogenetically related A. hydrophila and A. dhakensis, as well as between the phylogenetically distant A. hydrophila and A. salmonicida. Genes such as aerA, fla, laf, and lip were found across nearly all species. For instance, aerA was present in over 60% of A. caviae, A. dhakensis, and A. veronii, and in over 95% of A. hydrophila and A. salmonicida genomes. Similarly, fla was widely distributed, occurring in 70–90% of genomes across most species. Consistent with previous studies61,62, hlyA and ast were more species-specific, with hlyA almost exclusively found in A. dhakensis (99.5%) and A. hydrophila (100%), and ast present in 95.2% of A. hydrophila genomes but nearly absent elsewhere. A. caviae (n=608) and A. salmonicida (n=176) carried the highest proportions of key virulence genes; all A. caviae genomes encoded alt and lip, with nearly all (99.7%) possessing ela, while A. salmonicida uniformly carried aerA, alt, ela, lip, and act. Genes ela and lip were also prevalent in A. dhakensis (100%), A. hydrophila (100%) and A. sanarellii (93.8%), while A. enteropelogenes and A. jandaei lacked several virulence genes (Fig. 1; summarised in Fig. 3a).

The Type III secretion system genes (ascF-G) were found more often in A. salmonicida (80%, 141/176), A. veronii (51%, 288/566), A. dhakensis (42%, 78/186), and A. hydrophila (17%, 21/124). Gene hlyA was present in all isolates of A. salmonicida, A. dhakensis, and A. hydrophila, but absent in A. jandaei, A. enteropelogenes, A. sanarellii, A. caviae, and A. veronii, representing a potential molecular marker for differentiating between species (Fig. 2b). Similarly, the enterotoxin-encoding ast gene was found only in A. hydrophila (95%, 118/124) and at low levels in A. dhakensis (6%, 11/186), highlighting its potential for species-specific identification. We noted, a high prevalence of the combination of the cytotoxic enterotoxin (act) and cytotoxic heat-labile enterotoxin (alt) in A. salmonicida (alt=100%, 176/176; act=99%, 175/176), A. hydrophila (alt=100%, 124/124; act=72%, 90/124), and A. dhakensis (alt=94%, 176/186; act=29%, 55/186). We also observed variation in the distribution of laf, suggesting that this gene had been gained or lost multiple times across the phylogeny (Fig. 1); for example, it was present in A. salmonicida ST-2 but absent from ST-1799.

Lastly, we investigated whether virulence gene presence varied within species depending on sample origin. For this, we focused on South Asia (Bangladesh, India and Pakistan), since here we had 1438 Aeromonas spp. genomes confirmed to originate from clinical or environmental samples. As shown in Supplementary Fig. 5a and 5b, no obvious association was observed between species, sample origin, and virulence gene presence. This was also true when looking at the phylogeny of A. caviae taken from South Asia (n=577), where we had almost equal representation from clinical (n=286; of which 153 were from individuals with moderate to severe diarrhoea, while 133 were asymptomatic) and environmental (n=291) Aeromonas genomes (Supplementary Fig. 6). Of the 11 BAPS clusters, 10 contained both clinical and environmental Aeromonas genomes; the only exception was cluster-4.10, which included just five isolates in total, all from environmental sources. Sublineages containing a mixture of both clinical and environmental isolates were also identified in the other abundant Aeromonas species (A. veronii, A. dhakensis, A. enteropelogenes, and A. hydrophila; Supplementary Table 3).

To understand if there were differences in AMR profiles across or within species, we inferred the genotypic AMR profiles for all the genomes included in this study. Overall, we found 162 discrete AMR genes, conferring resistance to 16 different antimicrobial classes, including β-lactams, sulphonamides, tetracycline, and aminoglycosides (Fig. 3b, summarised in Supplementary Table 4). Of the 162 AMR genes identified, those conferring resistance to β-lactam antibiotics were the most common, representing 40.74% (66/162) of the genes.

Within this group, Class D oxacillinases were the most common, present in 94.5% of genomes, with genes such as blaOXA-427 (44.6%, 826/1853) and blaOXA-12 (32.7%, 605/1853) being particularly abundant. This was followed by Class C cephalosporinases, which were frequently found in 71.3% of genomes, including blaMOX-6 (18.9%, 350/1853) and blaFOX-2 (9.6%, 177/1853). Class B metallo-β-lactamases were detected in 55.4% of genomes, dominated by blaCEPH-A3 (22.1%, 410/1853) and blacphA4 (9.6%, 177/1853). In contrast, Class A β-lactamases were less common, identified in 5.07% of genomes, with genes such as blaKPC-1 (0.8%, 15/1853) and blaRSA-1 (1.1%, 20/1853) among the more notable types (Supplementary Table 5). Of note, mobile colistin resistance genes (mcr) were detected in several species. These included A. jandaei (100%, 25/25), all of which carried mcr-7.1; A. piscicola (100%, 1/1) harbouring mcr-3 and mcr-3.8; A. salmonicida (59.1%, 104/176) carrying both mcr-3 and mcr-3.12; A. media (28.6%, 2/7) with mcr-3 and mcr-3.6; A. sobria (16.7%, 1/6) harbouring mcr-3 and mcr-3.12; A. veronii (1.94%, 11/566) with mcr-3, mcr-3.12, mcr-3.3, mcr-3.6, and mcr-3.8; A. caviae (1.32%, 8/608) carrying mcr-3, mcr-3.12, mcr-3.3, and mcr-3.6; and A. hydrophila (0.81%, 1/124), which harboured mcr-5 (Supplementary Table 4; see Supplementary Data 1). Only three genomes in our collection lacked any known AMR gene.

Looking more closely at the subset of 1438 South Asian genomes, we found that genotypic AMR patterns varied significantly across India, Bangladesh, and Pakistan (Fig. 4a). Aeromonas isolates from Pakistan exhibited the highest prevalence of resistance to β-lactams (cephamycin: 65.8%, 290/441), sulphonamides (24.3%, 107/441), diaminopyrimidines (22.4%, 99/441), aminoglycosides (22.2%, 98/441), tetracyclines (18.4%, 81/441), phenicols (16.3%, 72/441), and fluoroquinolones (12.9%, 57/441). Indian isolates showed comparatively lower prevalence of resistance gene presence amongst the three South Asian countries, with the highest resistance observed for carbapenems (69.2%, 553/799); macrolides (7.8%, 62/799); rifamycins (3.3%, 26/799); and oxazolidinones, pleuromutilins, and streptogramins (4.4%, 35/799). Aeromonas from Bangladesh had the lowest prevalence of resistance, with no genes conferring resistance to lincosamides, oxazolidinones, pleuromutilins, streptogramins, and nucleosides. Notably, all genomes across the three countries were predicted to carry cephalosporin resistance genes. Among these, 4.5% of the isolates (65/1438; 14 clinical and 51 environmental) carried genes encoding third-generation extended-spectrum β-lactamases (ESBLs), such as blaCTX-M, blaGES, blaPER, blaTEM, and blaVEB.

To further refine our understanding beyond drug class-level resistance, we assessed the genetic differences in AMR gene profile between clinical and environmental isolates among the 1438 genomes using PCA and NMDS. This revealed considerable overlap between both groups within each species (Fig. 4c and Supplementary Fig. 7a). Notably, species like A. caviae and A. veronii formed distinct clusters, likely due to differences in AMR gene diversity (Supplementary Fig. 7b and Supplementary Table 6). Environmental isolates generally exhibited a greater diversity in AMR genes carriage than the clinical isolates across Aeromonas species (Supplementary Table 7), except in A. caviae, where clinical isolates contained more unique determinants (18 vs. 7) but also had largest overlap (48 shared genes; Supplementary Fig. 8). Notably, A. veronii and A. dhakensis environmental isolates carried extensive AMR repertoires (45 and 22 genes, respectively) with several shared with the clinical strains, suggesting potential environmental reservoirs (Supplementary Table 7). Shared genes across species commonly included tet(A), tet(E), blaOXA and blaMOX variants, sul1/sul2, and qnr variants, indicating frequent exchange between strains occurring in the aquatic and clinical settings, while certain ESBLs and carbapenemases (e.g., blaCTX-M-15, blaVIM-6) remained confined to the clinical isolates included in this study. Furthermore, certain AMR genes exhibited notable differences in prevalence between the two groups (Fig. 4d); clinical isolates had a higher proportion of blaOXA-427 (64.7%, 288/445), blaMOX-6 (35.7%, 159/445), and sul1 (17.3%, 77/445), compared to environmentally derived isolate genomes. In contrast, environmental isolates showed higher proportions of blaOXA-12 (43.9%, 435/990), blaOXA-724 (18.3%, 181/990), blaCEPHA3 (28.7%, 284/990) and blacphA4 (12.7%, 126/990) (Fig. 4d). Of note, one Indian environmentally derived A. sanarellii isolate was predicted to be multidrug resistant by possessing genes to seven drug classes including: aminoglycoside, β-lactam, diaminopyrimidine, macrolide, phenicol, rifamycin, and sulphonamide (see Supplementary Data 1).

Methods
We used a collection of 996 novel Aeromonas isolates. Of these, 198 had been originally cultured and identified as V. cholerae using conventional microbiological and biochemical techniques90 before they were sequenced and classified genomically as Aeromonas. These included 133 isolates from a study entitled the Cholera-Hospital-Based-Intervention-for-7-Days (ChoBI7) programme in Dhaka city49 and 65 isolates from an environmental surveillance study for V. cholerae conducted in Bangladesh. These 198 genomes were obtained from diverse sources, including drinking water, plankton, zooplankton, sediment, and human stool or rectal swabs, collected between 2004 and 2016. The remaining 798 genomes were obtained from 408 water samples collected for isolation of V. cholerae from ponds, rivers, and canals across six different states—Delhi, Haryana, Himachal Pradesh, Punjab, Uttar Pradesh and Uttarakhand, as well as Chandigarh in Northern India—between 2020 and 2023. The collected water samples were first enriched in alkaline peptone water and Selenite F broth. They were then sub-cultured onto various media. Blood agar and Thiosulphate Citrate Bile Salt Sucrose (TCBS; Difco Laboratories, Detroit, Michigan, USA) agar were used for the isolation of V. cholerae and Aeromonas spp. MacConkey agar (HiMedia Laboratories Private Limited, Mumbai, India) and Xylose Lysine Deoxycholate (XLD; Difco Laboratories, Detroit, Michigan, USA) were used for isolating other Enterobacteriaceae. The presumptively selected colonies from the media listed above were identified by matrix-assisted laser desorption ionisation-time of flight mass spectrometry (MALDI-TOF)91, and isolates that were identified as Aeromonas spp. underwent DNA extraction and whole genome sequencing as described below. Between 1 and 9 colonies were retained from each Aeromonas-positive water sample in India for downstream processing, including DNA extraction and whole genome sequencing.

DNA was extracted from bacterial isolates using the Wizard® Genomic DNA Purification Kit (Promega, Madison, WI, USA) and Qiagen DNeasy Blood & Tissue Kits (Qiagen, Hilden, Germany). Illumina sequencing libraries were prepared using 0.5µg DNA and sheared using an S2 ultrasonicator (Covaris), before adaptor ligation and dual index barcoding (NEBNext UltraII, New England Biolabs)92. Libraries were sequenced on the HiSeq X10 and NovaSeq 6000 platform (Illumina, San Diego, CA, USA) at the Wellcome Sanger Institute (Cambridge, UK), generating 150 bp paired-end reads with a mean depth of coverage of 68× per genome (range: 21–118×, excluding outliers).

A total of 1853 genomes were analysed, including 996 genomes sequenced in this study. These genomes were contextualised using 416 Aeromonas genomes (including one genome from India) downloaded from a library of 661,405 previously assembled bacterial genomes held in the European Nucleotide Archive (ENA)50, along with 441 genomes from a study on children with moderate to severe diarrhoea and corresponding controls in Pakistan46.

We performed initial quality control of sequencing reads using FastQC (v0.11.4.)93 and MultiQC (v1.17)94. To further screen for contamination in the raw sequencing reads, we used Kraken2 (v2.2)95 and Bracken (v2.6.2)96 (Standard RefSeq database containing archaea, bacteria, viruses, plasmids and human, downloaded from https://benlangmead.github.io/aws-indexes/k2 on 9th October 2023). No trimming of raw reads was performed prior to de novo assembly using SPAdes (v3.9.0)97 to obtain contiguous scaffolded assemblies. Additional quality control was performed on the de novo assemblies using Quast (v5.0.2)98 and CheckM (v1.1.2)99 to calculate the genome contamination and completeness. Assemblies were filtered to retain high-quality assemblies based on genome length (4.5-5.2 Mb), genome contamination (<5%) and genome completeness (>95%). Overall, 5% (97/1950) of the genomes were excluded from the study based on these criteria, including six genomes from the study in Pakistan.

For species-level classification, we employed GTDB-Tk (v.2.4.1)100 with the Genome Taxonomy Database (GTDB; R226, released on 16 April 2025), see Supplementary Data 2. GTDB-Tk integrates marker gene analysis, phylogenetic placement, and average nucleotide identity (ANI) to provide robust and consistent taxonomic assignments. This was followed by assessment of genome-relatedness indices, where we used FastANI (v1.3)101 for the pairwise comparison of the Aeromonas genomes. For a subset of genomes (n=49), digital DNA–DNA hybridisation (dDDH) assays were also performed to further refine species-level resolution (Supplementary Data 3). The FastANI threshold was set at 90% based on pairwise sequence mapping101, and heatmaps based on estimated ANI for all Aeromonas genomes were created using the Complexheatmap package in R (v4.4.3)102.

AMR genes were identified from the assembled contigs using ABRicate (v1.0.1, available at https://github.com/tseemann/abricate) to query the Comprehensive Antimicrobial Resistance Database (CARD, v3.2.4, last updated in August 2022)103, where the CARD identification threshold was set at ≥80% and <80% nucleotide identity, both with ≥80% coverage (Supplementary Data 1). For the identification of the putative virulence genes, we initially used ABRicate with the Virulence Finder Database104, but this did not yield any hits. We subsequently used in_silico_pcr (v1.0.0, available at https://github.com/sanger-pathogens/sh16_scripts/blob/master/legacy/in_silico_pcr.py) with previously described primer sequences (see Supplementary Table 8)105,106 to detect virulence genes. To explore the patterns of AMR and virulence gene distribution in relation to sample origin (clinical vs. environmental) and Aeromonas species, we conducted a Principal Component Analysis (PCA) using the prcomp function from the stats package in R (v4.4.3)102, based on presence/absence matrices. The first two components, PCA1 and PCA2, were used. Furthermore, we applied Non-metric Multidimensional Scaling (NMDS) with the metaMDS function from the vegan package in R (v4.4.3)102, to investigate whether the AMR gene diversity of sample clusters based on the species or the source of isolation (clinical and environmental), using Bray-Curtis dissimilarity to achieve higher resolution. The two-dimensional ordinations from both PCA and NMDS were visualised using the ggplot2 package.

Draft Aeromonas genomes were annotated using PROKKA (v1.4.5)107 with default parameters. We used Panaroo (v1.3.0)108 to infer the pangenome of 1853 genomes and used the resulting core gene alignment (comprising 1969 core genes) for phylogenetic analysis. We extracted 1,282,613 variable sites using SNP-sites (v2.5.1)109, and the resulting multiple sequence alignment file was used to infer a phylogeny with IQ-TREE (v2.1)110, using a generalised time reversible substitution model with a FreeRate model of heterogeneity111, incorporating invariant sites using the -fconst parameter, and with 1000 UltraFast Bootstraps112. Phylogenetic trees and associated metadata were visualised in iTOL (v.5.6)113. Similarly, a separate phylogenetic tree was constructed for genomes from India, Bangladesh, and Pakistan (n=1438), consisting of 2067 core genes and 1,197,118 variable sites.

We used the MLST tool (v2.9, available at https://github.com/tseemann/mlst) to infer MLST against the Aeromonas scheme developed by Martino and colleagues51 and hosted at https://pubmlst.org/organisms/aeromonas-spp114. This MLST scheme uses six housekeeping genes (gyrB, groL, gltA, metG, ppsA, and recA) to subtype Aeromonas below the species level into Sequence Types (STs). Novel alleles and STs were submitted to the Aeromonas MLST database via pubmlst.org for naming.

Given the broad distribution of genomes across diverse MLST types among species, we also employed Bayesian Analysis of Population Structure (BAPS) using fastBAPS (v1.0.6)115 to define higher-level groupings.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript investigates the genomic characteristics of Aeromonas spp. using a broad collection of isolates from South Asia. Although the dataset is extensive, the study contains several methodological and interpretative deficiencies that significantly limit the validity of its conclusions and should be thoroughly addressed before publication.\n\nThe Introduction section should not solely focus on human infections. Aeromonas spp. are also capable of causing infections in terrestrial and aquatic animals; this aspect should be elaborated. Moreover, the zoonotic potential of these species should be emphasized.\n\nA dedicated paragraph discussing the diagnostic challenges associated with Aeromonas identification and methods for achieving definitive diagnosis should be included in the Introduction.\n\nAdditional details should be provided for the isolates mentioned on lines 142–143, including the time frame of isolation and specific sources.\n\nOn line 442, the number of isolates subjected to genome analysis should be explicitly stated.\n\nAlthough the WGS platform (HiSeq X10) is mentioned, essential methodological details such as library preparation procedures, targeted sequencing depth, and genome coverage are lacking and should be added.\n\nThe methodology does not clearly indicate whether raw reads were subjected to trimming or quality filtering post-NGS. This information should be clarified.\n\nIn the genome-based identification section, it would be more appropriate to use DNA-DNA hybridization (dDDH) values instead of Kraken and Bracken. dDDH, especially when combined with ANI values, provides a more reliable taxonomic resolution.\n\nThe database used for species identification is from 2023. However, even within the current year, several new Aeromonas species have been described. The study should ensure the use of the most recent taxonomic references.\n\nIn the section on prediction of AMR and virulence genes, the threshold values (identity and coverage) should be clearly stated. Additionally, the dates of the last database updates should be provided.\n\nThe PubMLST link cited on line 497 is incorrect. It should be updated to the correct URL: \"https://pubmlst.org/organisms/aeromonas-spp\". The section titled \"Multi-locus Sequence Typing (MLST) and Subtyping of Aeromonas isolates\" should be revised accordingly.\n\nThe Results section lacks proper evaluation due to incorrect species assignments, which compromises the interpretation of findings.", "reviewer_2": "The manuscript presents a comprehensive view of the Aeromonas genus, including 1,853 genomes, and a detailed comparison of clinical and environmental strains from South Asia, with 996 newly sequenced genomes from Bangladesh and India. The authors have identified a high diversity of Aeromonas without a clear separation between clinical and environmental isolates. They also observed a high incidence of AMR genes across all isolates, including resistance to front-line and last-line antibiotics. Additionally, the manuscript highlights the frequent misidentification of Aeromonas as Vibrio cholerae, which is a key issue in cholera-endemic regions where both genera coexist and are associated with diarrheal disease. The findings are of interest to readers, but the representativeness of the dataset and the depth of the analysis are insufficient. Below are my comments:\n\n1.The title should be revised to explicitly state the geographic focus of the study, as the dataset primarily includes strains from South Asia. This is important for clarity and to avoid misleading readers about the scope of the research. Additionally, in Figure 1, the species Aeromonas salmonicida appears to be represented by nearly identical clones, which raises questions about the representativeness of the sampled strains. This might be a limitation in the study's ability to reflect the true genetic diversity of the genus.\n\n2.The authors only searched the ENR database for genome data, but the NCBI database currently contains over 3,000 Aeromonas genomes. Including these additional genomes would provide a more comprehensive view of the genus and enhance the understanding of the public health significance of Aeromonas in terms of AMR dissemination and global transmission dynamics.\n\n3.In the introduction section, it is recommended to appropriately introduce the current status of antibiotic resistance in Aeromonas species to enhance the reader's understanding. For example, existing studies have shown that Aeromonas is a natural host for MCR-3 and MCR-7, [10.1016/j.drup.2024.101088] [10.1016/j.drup.2023.101006].\n\n4.The manuscript's genome analysis lacks sufficient depth and accuracy, which limits the ability to draw clear conclusions. For example, the identification threshold for resistance genes was set at 80% nucleotide identity, but many studies have shown that resistance genes such as mcr-7 in Aeromonas often have lower nucleotide identity. Additionally, many resistance genes are naturally carried by Aeromonas. Whether these genes can evolve into acquired drug resistance genes and spread remains to be further analyzed.\n\n5.It is recommended that the authors supplement the manuscript with data on the phenotypic drug resistance of the strains. This would help clarify the contribution of the identified resistance genes to the overall drug resistance profile of the isolates.", "reviewer_3": "Synopsis\nThis study conducted environmental surveillance of Aeromonas spp. in India and Bangladesh and performed WGS on 996 isolates. They then combined these newly sequenced data with publicly available genomes of 441 clinical Aeromonas isolates from Pakistan, as well as a global collection of 416 Aeromonas isolates from different sources. A set of bioinformatic tools, including Kraken 2, Bracken, GTDB-Tk and FastANI were used to determine Aeromonas species. Standard genomic and phylogenetic tools were employed to define the population genetic structure and characterize the distribution of virulence and AMR genes by species and sample type.\n\nThey found a highly diverse Aeromonas genus with 28 species and 905 novel sequence types (STs), comprising 72.5% of total genomes. Clinical and environmental isolates clustered together across the phylogenetic tree. A high prevalence of AMR genes was found across isolates, including predicted resistance to last-line antibiotics. They further highlighted a frequent misidentification of Aeromonas as Vibrio cholera in cholera endemic settings. They concluded that Aeromonas is an important environmental AMR reservoir, encompassing multi-species capable of causing human infections.\n\nStrengths:\nThere has been a paucity of genomic data on Aeromonas, particularly in South Asian countries, where little is known about its environmental circulation and potential for spillover into human disease, such as childhood diarrhea. This study presents large-scale genomic surveillance of environmental isolates, offering new insights into the high species and genetic diversity of Aeromonas, its role as a major AMR reservoir, including resistance to last-resort antibiotics, and the strong genetic links between environmental and diarrhea-causing human isolates. Within this genomic framework, the authors have contributed to the expansion of the Aeromonas MLST database, with 905 novel sequence types now defined, providing a valuable resource for the broader research community. Notably, the high co-occurrence of Aeromonas and Vibrio cholerae in environmental samples is an important finding, raising concerns about the potential for overestimation of cholera prevalence in environmental surveillance, particularly in endemic settings.\n\nLimitations:\nBecause the authors compared species, sequence types, AMR, and virulence gene distributions between clinical isolates (sourced from Pakistan) and environmental isolates (from Bangladesh and Pakistan), and given that the sampling strategies differed between India (multiple colonies per sample) and Bangladesh (single colony per sample), the manuscript should include a paragraph discussing potential biases introduced by these differences in sampling approach, geographical coverage, and isolate sources. Additionally, the study would be strengthened by providing more detailed metadata regarding the environmental water sources; specifically, whether the samples were from drinking water, rivers, lakes, or streams. Information such as sampling locations, presence of antibiotic residues, linkage to human activity, and whether repeat sampling at the same sites occurred (if available) should be included.\n\nThroughout the manuscript, the authors use the term “clinical and environmental isolates.” It appears that this refers only to clinical isolates from Pakistan and environmental isolates from India and Bangladesh, excluding the global isolate collection. This definition should be made explicit to avoid confusion. Moreover, numerical values should be provided alongside all percentages mentioned in the manuscript to enhance clarity and interpretability.\n\nSpecific comments and questions:\n• Line 162: Please provide more details on the global isolates, including their sources (e.g., environment, human disease, animal disease, human carriage), whether they were associated with outbreaks, and their years of isolation.\n• Line 181: Consider rewriting and splitting this sentence for clarity.\n• Line 187: Provide the number of samples containing more than one species.\n• Line 188: Without associated metadata, interpreting country-specific species distribution from the global dataset is difficult. Please acknowledge this limitation or provide metadata.\n• Line 199: Do the authors have co-occurrence data for the Bangladesh site?\n• Line 207: Add a separate paragraph that describes in detail the phylogenetic relatedness of environmental and clinical isolates, particularly for species and specific sequence types where clusters are evident.\n• Line 229: Clarify the phrase \"human health.\" Also, consider to revise the header to specify \"Virulence Distribution.\"\n• Line 231: Add citations for each virulence factor mentioned.\n• Line 235: Can the authors explain why ABRicate with VFDB did not detect any hits? Did they attempt alternative tools like SRST2?\n• Line 239: The color coding in Figure 2a does not match the legend, please revise.\n• Line 240: Split the sentence into two and add more detail about the “virulence profile.”\n• Line 243: Specify what is meant by “others.”\n• Line 246: Clarify what is meant by “highest proportion of some virulence genes.” Which genes and which species?\n• Line 250: Add percentages where relevant.\n• Lines 251–253: Provide the proportions of hlyA and ast across species. If these genes are proposed as species markers, this information is essential.\n• Line 254: This sentence should be rewritten as two sentences for clarity.\n• Lines 260–265: The association between species and sample origin appears to have already been discussed in line 178 —please check for redundancy.\n• Line 266: This content would be better placed under the “Genomic Taxonomy” section header.\n• Line 274: Provide percentages for each β-lactam antibiotic subgroup.\n• Line 276: Provide the proportion of mcr genes found in each of the eight species. Were these isolates phenotypically tested for colistin resistance?\n• Line 282: It seems that the authors are comparing AMR gene distributions across countries, not phenotypic resistance,please clarify. Note that the correlation between AMR genotype and phenotype is not well-studied in these organisms.\n• Lines 291–296: These sentences would be better placed in the Discussion and expanded upon.\n• Line 298: Consider including comparisons of AMR gene profiles between clinical and environmental isolates for the predominant species.\n• Line 324: Consider adding more references (e.g., PMID: 38739115; PMID: 37808293).\n• Line 326: According to Supplementary Figure 1B, A. dhakensis may not be overpresented in Bangladesh, please clarify.\n• Lines 327–331: Metadata for global isolates are needed to understand their epidemiological context.\n• Line 346: Is this line referring to diarrheal disease? In PMID: 38739115, A. dhakensis appear to be more prevalent than other Aeromonas species in human bloodstream infections.\n• Line 353: Expand on the interpretation of the presented genetic data.\n• Line 354: Please provide further clarification, particularly if this information is intended to inform the development of future assays for Aeromonas species classification\n• Line 364: Confirm whether a substantial number of ESBL genes were identified.\n• Line 365: Replace with “AMR gene profile.”\n• Lines 365–374: Differences in AMR and virulence gene distribution may reflect the fact that all clinical isolates were from Pakistan, while environmental isolates were more diverse and came from two countries. This should be discussed.\n• Line 387: Please provide more explanation on why isolates from Pakistan harbored a higher number of resistance genes. Could this be attributed to their clinical origin, where selection pressure from antibiotic use is typically higher? Alternatively, could the difference be influenced by the absence of environmental isolates from Pakistan in the study, limiting a direct comparison across sample types?\n• Line 427: Please provide additional metadata for the isolates from Bangladesh.\n• Line 501: This sentence is difficult to understand,please rewrite for clarity.\n• Line 807: Consider revising to: “Each dot is coloured based on the gene abundance per species (see scale), and the size of the dots corresponds to the number of genomes harboring the gene.”"}, "author_response": "TITLE: Aeromonas in South Asia: Genomic Insights into an Environmental Pathogen and\nReservoir of Antimicrobial Resistance.\nNature Communications manuscript NCOMMS-25-34237\nWe are sincerely thankful to the reviewers for their thoughtful and constructive comments. Their\nperspectives have helped us refine the manuscript in meaningful ways. Below, we provide detailed\nresponses to each point raised, and we have revised the manuscript accordingly.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author)\nThis manuscript investigates the genomic characteristics of Aeromonas spp. using a broad collection of\nisolates from South Asia. Although the dataset is extensive, the study contains several methodological\nand interpretative deficiencies that significantly limit the validity of its conclusions and should be\nthoroughly addressed before publication.\n1. The Introduction section should not solely focus on human infections. Aeromonas spp. are also capable\nof causing infections in terrestrial and aquatic animals; this aspect should be elaborated. Moreover, the\nzoonotic potential of these species should be emphasized.\nDone. Please see lines 109-114 for the updated text.\n2. A dedicated paragraph discussing the diagnostic challenges associated with Aeromonas identification\nand methods for achieving definitive diagnosis should be included in the Introduction.\nDone. Please see lines 132-141.\n3. Additional details should be provided for the isolates mentioned on lines 142–143, including the time\nframe of isolation and specific sources.\nDone. Please see lines 160-162.\n4. On line 442, the number of isolates subjected to genome analysis should be explicitly stated.\nDone. Please see line 513.\n5. Although the WGS platform (HiSeq X10) is mentioned, essential methodological details such as library\npreparation procedures, targeted sequencing depth, and genome coverage are lacking and should be\nadded.\nDone. Please see lines 508-511.\n6. The methodology does not clearly indicate whether raw reads were subjected to trimming or quality\nfiltering post-NGS. This information should be clarified.\nTo clarify, no trimming of raw reads was performed prior to de novo assembly. In addition, CheckM\nwas also used for quality assurance to exclude incomplete or contaminated assemblies. This is all now\ndetailed in the Methods lines 520-524.\n7. In the genome-based identification section, it would be more appropriate to use DNA-DNA\nhybridization (dDDH) values instead of Kraken and Bracken. dDDH, especially when combined with\nANI values, provides a more reliable taxonomic resolution.\nWe apologise; this was wrongly stated in the Methods (previously line 462). Kraken/Bracken were used\nonly to screen for contamination, not for species delineation (see lines 521–523).\nSpecies identification was performed using GTDB-Tk (previously line 464, now lines 532–535), which\nprovides taxonomic assignments based on the Genome Taxonomy Database (GTDB), a widely\nrecognized reference standard for bacterial and archaeal genomes (Chaumeil et al., 2020, PMID:\n31730192). GTDB-Tk integrates average nucleotide identity (ANI), marker gene-based classification,\nand phylogenetic placement to provide a comprehensive species-level taxonomic assessment.\nTo further support species assignments, we applied FastANI, which provides ANI estimates comparable\nto dDDH (~70% dDDH = ~95% ANI; Goris et al., 2007, PMID: 17220447). Given this concordance\nand to avoid using multiple overlapping tools, we focused on GTDB-Tk with FastANI for efficient and\nconsistent species assignment in the manuscript.\n8. The database used for species identification is from 2023. However, even within the current year,\nseveral new Aeromonas species have been described. The study should ensure the use of the most recent\ntaxonomic references.\nIn response to your question, we repeated this analysis with the latest GTDB release (R226, 16 April\n2025; see Supplementary Data-2), which did not change the results for species diversity or composition.\n9. In the section on prediction of AMR and virulence genes, the threshold values (identity and coverage)\nshould be clearly stated. Additionally, the dates of the last database updates should be provided.\nDone. Please see lines 543 and 544.\n10. The PubMLST link cited on line 497 is incorrect. It should be updated to the correct URL:\n\"https://pubmlst.org/organisms/aeromonas-spp [pubmlst.org]\". The section titled \"Multi-locus\nSequence Typing (MLST) and Subtyping of Aeromonas isolates\" should be revised accordingly.\nDone. Please see lines 242 and 571.\n11. The Results section lacks proper evaluation due to incorrect species assignments, which compromises\nthe interpretation of findings.\nPlease see the responses to points 7 and 8 above. Given this we would argue that the species\nassignments are robust and therefore so too is the interpretations of the data.\nReviewer #2 (Remarks to the Author):\nThe manuscript presents a comprehensive view of the Aeromonas genus, including 1,853 genomes, and\na detailed comparison of clinical and environmental strains from South Asia, with 996 newly sequenced\ngenomes from Bangladesh and India. The authors have identified a high diversity of Aeromonas without\na clear separation between clinical and environmental isolates. They also observed a high incidence of\nAMR genes across all isolates, including resistance to front-line and last-line antibiotics. Additionally,\nthe manuscript highlights the frequent misidentification of Aeromonas as Vibrio cholerae, which is a\nkey issue in cholera-endemic regions where both genera coexist and are associated with diarrheal\ndisease. The findings are of interest to readers, but the representativeness of the dataset and the depth\nof the analysis are insufficient. Below are my comments:\n1. The title should be revised to explicitly state the geographic focus of the study, as the dataset primarily\nincludes strains from South Asia. This is important for clarity and to avoid misleading readers about the\nscope of the research.\nThank you for your comment this was extremely helpful. As suggested, we have revised the title.\n2. Additionally, in Figure 1, the species Aeromonas salmonicida appears to be represented by nearly\nidentical clones, which raises questions about the representativeness of the sampled strains. This might\nbe a limitation in the study's ability to reflect the true genetic diversity of the genus.\nThe observation from Figure 1 that Aeromonas salmonicida appears to be represented by nearly\nidentical clones, we would like to highlight that among the 176 A. salmonicida genomes, the majority\noriginate from Denmark (100 genomes, 1980–2011) and China (60 genomes, 2012–2015), covering a\nbroad temporal range. The remaining genomes come from other countries. A comprehensive sampling\noverview has been added to the revised manuscript to further clarify the dataset’s composition and its\nimplications for genetic diversity (see Supplementary Table 1).\n3. The authors only searched the ENR database for genome data, but the NCBI database currently contains\nover 3,000 Aeromonas genomes. Including these additional genomes would provide a more\ncomprehensive view of the genus and enhance the understanding of the public health significance of\nAeromonas in terms of AMR dissemination and global transmission dynamics.\nLikely because of the increasing significance of Aeromonas in human infection, this is a rapidly\nchanging situation with more genomes are appearing each month. To address your comment, we\nretrieved all 4,378 Aeromonas genomes from NCBI (as of August 18, 2025) and classified them using\nGTDB-Tk. After excluding 30 assemblies that could not be assigned at the species level, 4,344\nassemblies remained.\nOverall, the NCBI dataset contained 32 Aeromonas species, the majority of which were represented by\nfive species (A. veronii, A. caviae, A. hydrophila, A. salmonicida, and A. dhakensis), which together\naccounted for ~86% (3,722/4,344) of genomes. This species distribution is, consistent with our data\n(~89%, 1,660/1,853). The species present in the NCBI dataset but absent from our study (A. aquatica,\nA. finlandensis, A. lacus, A. lusitana, and A. molluscorum) accounted for only ~0.5% of the total NCBI\ngenomes (20/4,344).\nHowever, NCBI dataset has a strong geographic bias, dominated by genomes from China (1,436), and\nfrom East Asia (1,785 genomes; ~41%), please see Figure below. Only 531 of the 4,344 genomes\n(~12%) originate from South Asia, (most of which are unpublished sequences from Bangladesh). Hence\nin principle this extra data is not relevant to the geographic focus of this study and given the 531 have\nonly recently been submitted and are unpublished we are reluctant to include them in our publication,\nin fear of disadvantaging others. However, to address your question and to increase the species and\ngeographic contextual data, without losing our focus on South Asia we incorporated the curated dataset\npublished by Grace Blackwell et al. (PMID: PMC8577725) to our phylogenies, which included 28\nAeromonas species from 16 countries across five continents and diverse sources (see revised\nSupplementary Table 1).\nFigure: Overview of Aeromonas genomes available in GenBank (as of August 18, 2025).\n4. In the introduction section, it is recommended to appropriately introduce the current status of antibiotic\nresistance in Aeromonas species to enhance the reader's understanding. For example, existing studies\nhave shown that Aeromonas is a natural host for MCR-3 and MCR-7, [10.1016/j.drup.2024.101088]\n[10.1016/j.drup.2023.101006].\nThank you for the suggestion. We have revised the Introduction to include this information. Please see\nlines 143-151.\n5. The manuscript's genome analysis lacks sufficient depth and accuracy, which limits the ability to draw\nclear conclusions. For example, the identification threshold for resistance genes was set at 80%\nnucleotide identity, but many studies have shown that resistance genes such as mcr-7 in Aeromonas\noften have lower nucleotide identity. Additionally, many resistance genes are naturally carried by\nAeromonas. Whether these genes can evolve into acquired drug resistance genes and spread remains to\nbe further analyzed.\nThank you for your comment. We reanalysed the data using a more inclusive threshold of ≥30%\nnucleotide identity combined with ≥80% coverage. Importantly, below the 80% nucleotide identity\nthreshold, only the following MCR variants were identified: MCR-3, MCR-3.6, MCR-3.8, MCR-3.9,\nMCR-3.12, MCR-7.1, and MCR-9; no other AMR genes were detected. This information is presented\nin Supplementary Data-1, under a separate column named “AMR_Genes (Coverage ≥80%, Nucleotide\nIdentity <80%).”\n6. It is recommended that the authors supplement the manuscript with data on the phenotypic drug\nresistance of the strains. This would help clarify the contribution of the identified resistance genes to\nthe overall drug resistance profile of the isolates.\nHere, unfortunately phenotypic AST testing was not standard practise or standardised across the\nmultiple different labs collecting these isolates, especially for environmentally derived Aeromonas\nisolates. What is more, we did not have access to the live isolates, hence we could not test them\nconsistently either. Given this and to ensure consistency we chose to focus on the genetic differences,\nincluding in antimicrobial resistance gene profiles, between Aeromonas isolates from different sample\nsites and across species. Given the above and whilst these data provide important insights into the\nresistome and the associated risks to both aquatic environments and human health we would argue that\nadding this analysis is out of scope for this study. We have now clarified this limitation in the Discussion\nsection (Please see lines 448-449).\nReviewer #3 (Remarks to the Author)\nSynopsis\nThis study conducted environmental surveillance of Aeromonas spp. in India and Bangladesh and\nperformed WGS on 996 isolates. They then combined these newly sequenced data with publicly\navailable genomes of 441 clinical Aeromonas isolates from Pakistan, as well as a global collection of\nBracken, GTDB-Tk and FastANI were used to determine Aeromonas species. Standard genomic and\nphylogenetic tools were employed to define the population genetic structure and characterize the\ndistribution of virulence and AMR genes by species and sample type.\nThey found a highly diverse Aeromonas genus with 28 species and 905 novel sequence types (STs),\ncomprising 72.5% of total genomes. Clinical and environmental isolates clustered together across the\nphylogenetic tree. A high prevalence of AMR genes was found across isolates, including predicted\nresistance to last-line antibiotics. They further highlighted a frequent misidentification of Aeromonas\nas Vibrio cholera in cholera endemic settings. They concluded that Aeromonas is an important\nenvironmental AMR reservoir, encompassing multi-species capable of causing human infections.\nStrengths:\nThere has been a paucity of genomic data on Aeromonas, particularly in South Asian countries, where\nlittle is known about its environmental circulation and potential for spillover into human disease, such\nas childhood diarrhea. This study presents large-scale genomic surveillance of environmental isolates,\noffering new insights into the high species and genetic diversity of Aeromonas, its role as a major AMR\nreservoir, including resistance to last-resort antibiotics, and the strong genetic links between\nenvironmental and diarrhea-causing human isolates. Within this genomic framework, the authors have\ncontributed to the expansion of the Aeromonas MLST database, with 905 novel sequence types now\ndefined, providing a valuable resource for the broader research community. Notably, the high co-\noccurrence of Aeromonas and Vibrio cholerae in environmental samples is an important finding, raising\nconcerns about the potential for overestimation of cholera prevalence in environmental surveillance,\nparticularly in endemic settings.\nLimitations:\nBecause the authors compared species, sequence types, AMR, and virulence gene distributions between\nclinical isolates (sourced from Pakistan) and environmental isolates (from Bangladesh and Pakistan),\nand given that the sampling strategies differed between India (multiple colonies per sample) and\nBangladesh (single colony per sample), the manuscript should include a paragraph discussing potential\nbiases introduced by these differences in sampling approach, geographical coverage, and isolate\nsources. Additionally, the study would be strengthened by providing more detailed metadata regarding\nthe environmental water sources; specifically, whether the samples were from drinking water, rivers,\nlakes, or streams. Information such as sampling locations, presence of antibiotic residues, linkage to\nhuman activity, and whether repeat sampling at the same sites occurred (if available) should be included.\nThroughout the manuscript, the authors use the term “clinical and environmental isolates.” It appears\nthat this refers only to clinical isolates from Pakistan and environmental isolates from India and\nBangladesh, excluding the global isolate collection. This definition should be made explicit to avoid\nconfusion. Moreover, numerical values should be provided alongside all percentages mentioned in the\nmanuscript to enhance clarity and interpretability.\nSpecific comments and questions:\n1. Line 162: Please provide more details on the global isolates, including their sources (e.g., environment,\nhuman disease, animal disease, human carriage), whether they were associated with outbreaks, and their\nyears of isolation.\nThank you for your comment. We have included the available metadata for the global isolates in the\nrevised manuscript. Please see lines 181-183 and Supplementary Table 1.\n2. Line 181: Consider rewriting and splitting this sentence for clarity.\nDone. Please see lines 201-204.\n3. Line 187: Provide the number of samples containing more than one species.\nDone. Please see line 206 and updated Supplementary Fig. 3a.\n4. Line 188: Without associated metadata, interpreting country-specific species distribution from the\nglobal dataset is difficult. Please acknowledge this limitation or provide metadata.\nDone. To address this comment Metadata has been included in the revised manuscript (please see\nSupplementary Table 1).\n5. Line 199: Do the authors have co-occurrence data for the Bangladesh site?\nUnfortunately, this is not available.\n6. Line 207: Add a separate paragraph that describes in detail the phylogenetic relatedness of\nenvironmental and clinical isolates, particularly for species and specific sequence types where clusters\nare evident.\nAdded, please see lines 296-298 and revised Supplementary Table 3.\n7. Line 229: Clarify the phrase \"human health.\" Also, consider to revise the header to specify \"Virulence\nDistribution.\"\nDone. We have revised the header to explicitly reflect the focus on virulence gene distribution (please\nsee line 250).\n8. Line 231: Add citations for each virulence factor mentioned.\nAdded, please see lines 252-255.\n9. Line 235: Can the authors explain why ABRicate with VFDB did not detect any hits? Did they attempt\nalternative tools like SRST2?\nWe had not used SRST2 in our original analysis. However, to address your comment, we ran SRST2\nand similarly found no virulence gene matches. ABRicate with VFDB did not detect any virulence genes\nbecause the VFDB database lacks curated entries specific to Aeromonas spp., focusing mainly\non Enterobacteriaceae and other well-characterized pathogens. Thus, it is clear that the absence of hits\nis due to the database content rather than the tool used. To address this limitation, we conducted a\nthorough literature review to identify experimentally validated Aeromonas-associated virulence genes.\nWe then compiled a custom gene set and used these data as a bespoke database with which to identify\nall known virulence genes across our genomes using in-silico PCR.\n10. Line 239: The color coding in Figure 2a does not match the legend, please revise.\nDone.\n11. Line 240: Split the sentence into two and add more detail about the “virulence profile”.\nDone. Additional details about the “virulence profile” have been included (please see lines 260-269).\n12. Line 243: Specify what is meant by “others”.\nAdded, please see lines 267.\n13. Line 246: Clarify what is meant by “highest proportion of some virulence genes.” Which genes and\nwhich species?\nDone. Please see lines 269-273.\n14. Line 250: Add percentages where relevant.\nDone.\n15. Lines 251–253: Provide the proportions of hlyA and ast across species. If these genes are proposed as\nspecies markers, this information is essential.\nDone.\n16. Line 254: This sentence should be rewritten as two sentences for clarity.\nDone. Please see lines 281-285.\n17. Lines 260–265: The association between species and sample origin appears to have already been\ndiscussed in line 178—please check for redundancy.\nThank you, the difference is that lines 260–265 (now lines 287–294) refer to the association between\nvirulence gene presence and sample of origin.\n18. Line 266: This content would be better placed under the “Genomic Taxonomy” section header.\nThank you, we have revised the section for clarity.\n19. Line 274: Provide percentages for each β-lactam antibiotic subgroup.\nDone. Please see lines 306-312.\n20. Line 276: Provide the proportion of mcr genes found in each of the eight species. Were these isolates\nphenotypically tested for colistin resistance?\nWe have updated the section to include mcr gene distribution per species. To clarify these isolates were\nnot phenotypically tested. Please also see response to Reviewer 2, question 6.\n21. Line 282: It seems that the authors are comparing AMR gene distributions across countries, not\nphenotypic resistance, please clarify. Note that the correlation between AMR genotype and phenotype\nis not well-studied in these organisms.\nDone.\n22. Lines 291–296: These sentences would be better placed in the Discussion and expanded upon.\nDone. Please see lines 413-420.\n23. Line 298: Consider including comparisons of AMR gene profiles between clinical and environmental\nisolates for the predominant species.\nTo address this, we have provided detailed AMR gene information for clinical and environmental\nisolates by species (please see lines 335-348, Supplementary Figures 7 and 8, and Supplementary Table\n6 and 7).\n24. Line 324: Consider adding more references (e.g., PMID: 38739115; PMID: 37808293).\nAdded.\n25. Line 326: According to Supplementary Figure 1B, A. dhakensis may not be overpresented in\nBangladesh, please clarify.\nCorrected, thank you for pointing this out.\n26. Lines 327–331: Metadata for global isolates are needed to understand their epidemiological context.\nAdded, please see lines 373-378.\n27. Line 346: Is this line referring to diarrheal disease? In PMID: 38739115, A. dhakensis appear to be\nmore prevalent than other Aeromonas species in human bloodstream infections.\nThank you for your comment. Yes, in line 346 we are referring to diarrheal cases, not bloodstream\ninfections. We had revised the text for clarity.\n28. Line 353: Expand on the interpretation of the presented genetic data.\nDone, please see lines 397-403.\n29. Line 354: Please provide further clarification, particularly if this information is intended to inform the\ndevelopment of future assays for Aeromonas species classification.\nDone, please see line 400-408.\n30. Line 364: Confirm whether a substantial number of ESBL genes were identified.\nDone, please see lines 416-417.\n31. Line 365: Replace with “AMR gene profile.\nDone.\n32. Lines 365–374: Differences in AMR and virulence gene distribution may reflect the fact that all clinical\nisolates were from Pakistan, while environmental isolates were more diverse and came from two\ncountries. This should be discussed.\nThank you for your comment. We have revised the discussion to make it clear that clinical isolates were\nall from Pakistan whereas the environmental isolates were from India and Bangladesh and noted that\nthis sampling bias may have influenced the observed differences in gene distribution (see revised text,\nlines 441–443).\n33. Line 387: Please provide more explanation on why isolates from Pakistan harbored a higher number of\nresistance genes. Could this be attributed to their clinical origin, where selection pressure from antibiotic\nuse is typically higher? Alternatively, could the difference be influenced by the absence of\nenvironmental isolates from Pakistan in the study, limiting a direct comparison across sample types?\nAdded, please see lines 443-451.\n34. Line 427: Please provide additional metadata for the isolates from Bangladesh.\nAdded, please see lines 492-494.\n35. Line 501: This sentence is difficult to understand, please rewrite for clarity.\nThank you, we have rephrased this sentence (please see lines 574-575).\n36. Line 807: Consider revising to: “Each dot is coloured based on the gene abundance per species (see\nscale), and the size of the dots corresponds to the number of genomes harboring the gene.”\nDone. Thank you"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Dear Authors,\nThank you for the substantial and carefully executed revision of “Aeromonas in South Asia: Genomic Insights into an Environmental Pathogen and Reservoir of Antimicrobial Resistance.” I appreciate the scope of the study and the evident effort of the team, and I offer the following comments purely in a constructive spirit.\n\nMy remaining concern is the taxonomic identification of the Aeromonas isolates. Based on my research and experience with Aeromonas, average nucleotide identity (ANI) alone is not sufficient for species-level assignments, particularly in borderline cases. I strongly recommend that you combine ANI with digital DNA–DNA hybridization (dDDH), reporting both metrics and using accepted thresholds (e.g., ANI ≈95–96% and dDDH ≈70%) to reach stable conclusions. Please provide a table (preferably supplementary) listing, for each isolate, the closest type strain(s) together with ANI and dDDH values and the corresponding genome accessions, along with the software and parameters used (e.g., FastANI/pyani for ANI; GGDC/Formula 2 for dDDH).\n\nRelatedly, while GTDB/GTDB-Tk is valuable for genome-based taxonomy, it is not a nomenclatural authority and may include provisional or non-validated names and an incomplete set of type strains. In the current GTDB release you cite, I noted ~32 Aeromonas taxa irrespective of nomenclatural status, whereas the List of Prokaryotic names with Standing in Nomenclature (LPSN; https://lpsn.dsmz.de/genus/aeromonas\n) currently recognises 34 validly published Aeromonas species. Please make sure your reference panel includes the type strain genomes of all valid species listed by LPSN, not only those represented in GTDB. Augmenting the GTDB-Tk reference with missing valid Aeromonas type strains may change some assignments, and this possibility should be checked and, if relevant, discussed.\n\nAddressing these points will, I believe, make the taxonomic backbone of the manuscript more robust and increase confidence in the downstream ecological and AMR interpretations. Thank you again for your thorough work; I hope these suggestions are helpful.", "reviewer_2": "The authors have resolved all my issues.", "reviewer_3": "I thank the authors for their thorough revision of the manuscript and for addressing my previous comments and questions.\nThe data interpretation and overall clarity have been substantially improved. One minor comment remains: please revise the color scheme and legend of Supplementary Figure 5b to enhance clarity."}, "author_response": "TITLE: Aeromonas in South Asia: Genomic Insights into an Environmental Pathogen and Reservoir\nof Antimicrobial Resistance.\nNature Communications manuscript NCOMMS-25-34237A\nWe are grateful to the reviewers for their thoughtful and constructive comments. Below, we provide detailed\nresponses to each point raised and have revised the manuscript accordingly.\nREVIEWER COMMENTS IN REVERSE ORDER FOR CLARITY\nReviewer #3 (Remarks to the Author):\nI thank the authors for their thorough revision of the manuscript and for addressing my previous comments\nand questions. The data interpretation and overall clarity have been substantially improved. One minor\ncomment remains: please revise the color scheme and legend of Supplementary Figure 5b to enhance clarity.\nThank you for your comment, we have revised Supplementary Figure 5b for improved clarity.\nReviewer #2 (Remarks to the Author):\nThe authors have resolved all my issues.\nMany thanks for your suggestions on this manuscript.\nReviewer #1 (Remarks to the Author):\nDear Authors,\nThank you for the substantial and carefully executed revision of “Aeromonas in South Asia: Genomic\nInsights into an Environmental Pathogen and Reservoir of Antimicrobial Resistance.” I appreciate the scope\nof the study and the evident effort of the team, and I offer the following comments purely in a constructive\nspirit.\nThank you for your kind comments.\nMy remaining concern is the taxonomic identification of the Aeromonas isolates. Based on my research and\nexperience with Aeromonas, average nucleotide identity (ANI) alone is not sufficient for species-level\nassignments, particularly in borderline cases. I strongly recommend that you combine ANI with digital\nDNA–DNA hybridization (dDDH), reporting both metrics and using accepted thresholds (e.g., ANI ≈95–\n96% and dDDH ≈70%) to reach stable conclusions. Please provide a table (preferably supplementary) listing,\nfor each isolate, the closest type strain(s) together with ANI and dDDH values and the corresponding genome\naccessions, along with the software and parameters used (e.g., FastANI/pyani for ANI; GGDC/Formula 2\nfor dDDH).\nAs requested, the closest type strain(s) and corresponding ANI values for each isolate are included in\nSupplementary Data-2 as part of the GTDB-Tk output, which reports the nearest genome matches and\nassociated ANI values for all isolates. For dDDH values please see below.\nRelatedly, while GTDB/GTDB-Tk is valuable for genome-based taxonomy, it is not a nomenclatural\nauthority and may include provisional or non-validated names and an incomplete set of type strains. In the\ncurrent GTDB release you cite, I noted ~32 Aeromonas taxa irrespective of nomenclatural status, whereas\nthe List of Prokaryotic names with Standing in Nomenclature (LPSN; https://lpsn.dsmz.de/genus/aeromonas\n[lpsn.dsmz.de]) currently recognises 34 validly published Aeromonas species. Please make sure your\nreference panel includes the type strain genomes of all valid species listed by LPSN, not only those\nrepresented in GTDB. Augmenting the GTDB-Tk reference with missing valid Aeromonas type strains may\nchange some assignments, and this possibility should be checked and, if relevant, discussed.\nAddressing these points will, I believe, make the taxonomic backbone of the manuscript more robust and\nincrease confidence in the downstream ecological and AMR interpretations. Thank you again for your\nthorough work; I hope these suggestions are helpful.\nWe carefully considered the reviewer’s suggestion regarding the inclusion of digital DNA–DNA\nhybridization (dDDH) values. As recommended, we used the Type (Strain) Genome Server (TYGS/ GGDC\ndenoted TYGS from hence forth) linked from the LPSN web portal. Noting that TYGS allows only 50 genomes\nper submission and bootstrap-based confidence intervals are disabled when processing more than 20\ngenomes (see, GGDC FAQ here; hence not suitable for the 1,853 genomes analysed here), to address the\nreviewers comments directly we subsampled 50 phylogenetically representative genomes covering all\nspecies included in our study and assessed their dDDH values. Of the genomes analysed the species\nclassification was identical for 42 of the 50 genomes analysed on either TYGS or GTDB-Tk (Table-1):\nindicating a high concordance between GGDC and GTDB-Tk classifications.\nOf the remaining eight genomes for 6 of these, the closest matching species (A. veronii & A.\nallosaccharophila) from TYGS was the same as the one determined by GTDB-Tk. However, due to borderline\ndDDH scores (Aeromonas veronii CECT 4257; d4 = 67.7–69.9) they were either marked as unclassified or\npotential new species by TYGS (Table-1). Taking the borderline calls Hu et al. (2022,\nhttps://doi.org/10.3389/fmicb.2022.910277 ), suggest that when the dDDH value is close to, but below, the\n~70 % threshold, the taxonomic status of those strain should be confirmed using a combination of\nphenotypic, chemotaxonomic, and phylogenomic analyses. Given the two approaches above, and that all\nthese genomes had an ANI>95% to the closest reference by GTDB-Tk it appears that the overall species\nidentifications here were robust using the approaches described in our original analysis (consistent with\nGoris et al., 2007).\nThe 7th genome was only matched to the genus level by GTDB-Tk and assigned to ‘a potential new species’\nby TYGS. The 8th genome was discordant between the two approaches, likely because the A. rivipollensis\nreference genome is not available in TYGS (For a summary of the dual taxa calls see attached Table-1).\nIt is also important to note that many of those with dDDH (formula d4) values above 70% significantly\nmatched more than one Aeromonas species, despite the deep phylogenetic divisions: For example, strain\n‘23386_4#62’, originally classified as A. caviae, shows a dDDH (d4) value of 82.2% with A. caviae NCTC\n12244 and 81.3% with A. hydrophila subsp. anaerogenes CECT 4221 (Table-1). Given the limitations of\nscalability for TYGS and that our approach, combining GTDB-Tk and FastANI, provides a standardised,\nreproducible, and computationally efficient framework that has been widely used for bacterial taxonomy,\nwe respectfully argue not to include dDDH values into our analysis."}]

Decision letter