Five A. hydrophila isolates of varied host origin were used for reference purposes and these isolates originated from Thailand (isolate F2D20 from Rana rugulosa), Bangladesh (isolate T4 from Labeo rohita and isolate B2/12 from an unknown host), India (isolate VDS from Ictalurus punctatus), and USA (isolate AL09-71 from I. punctatus). A further reference isolate, Aeromonas veronii biovar sobria LMG 13067 (originating from an unspecified frog host in the USA), and 12 type strains of Aeromonas spp. were also included (Supplementary Table 1).

P. hypophthalmus exhibiting classical signs of MAS, including reddened fins and external and/or internal hemorrhage, were sampled from disease outbreaks in the Mekong Delta region (An Giang, Ben Tre, Can Tho, Dong Thap, Hậu Giang, Tien Giang, and Vinh Long provinces) and Dong Nai province during 2013–2019. For each sampling occasion, at least two representative field isolates were included that had been derived from separate farms located in the same geographic area. Additionally, nine farms from five regions were sampled more extensively between 2018 and 2019, where multiple fish, internal organs and isolates per culture plate were collected. After disinfecting the fish surface with 70% ethanol in water, the internal organs (liver, spleen, and head kidney) were dissected out aseptically, streaked on to Rimler-Shotts agar (RS; Himedia, India) or Aeromonas medium (supplemented with ampicillin; Thermo Fisher Scientific, UK) and incubated overnight at 28°C. Smooth, convex, round, and yellow (or dark green on Aeromonas medium) single colonies presumed to be Aeromonas spp. were sub-cultured in tryptone soya broth (TSB; Himedia) at 28°C for 16–20 h at 240 rpm. All isolates, samples and locations are listed in Supplementary Table 1. Cultures were cryopreserved at –80°C in TSB medium supplemented with 20% glycerol (Merck, USA).

Genomic DNA was extracted from the culture pellet of each isolate using the GeneJET Genomic DNA Purification kit (Thermo Scientific, USA). PCR was performed against the 16S rRNA gene using the primers of Trakhna et al. (2009), and this assay yields an amplicon of 103 bp in the presence of genomic DNA from A. hydrophila. This 16S rRNA PCR assay was chosen to identify A. hydrophila instead of the more traditional aerolysin gene PCR (Pollard et al., 1990), as initial evaluations found equivalent PCR outcomes with reference isolates within the Aeromonas genus whilst avoiding the presence of non-specific artifact products that appeared for a small proportion of aerolysin PCR-negative samples collected from the field (e.g., isolate TN120 from An Giang and isolates TN103 and TN117 from Dong Thap; data not shown). Each PCR contained 10 ng genomic DNA in a 10 μL total volume containing 1 × MyTaq HS PCR master mix (Meridian Bioscience, UK) and 0.2 μM of each primer. The amplification conditions were as follows: initial denaturation at 95°C for 3 min, followed by 30 cycles of 95°C for 15 s, 55°C for 15 s, and 72°C for 10 s. PCR amplicons were electrophoresed through 1% agarose gel (Meridian Bioscience, UK) containing 0.1 mg/mL ethidium bromide in 0.5 × TAE (20 mM Tris, 10 mM acetic acid, and 0.5 mM EDTA). The gel was observed under UV light using the InGenius gel imaging system (Syngene, UK) for the presence of the expected band.

Each field isolate was genotyped using a rep-PCR based on the single (GTG)₅ repetitive primer, as described previously (Bartie et al., 2012). Field isolates negative by 16S rRNA PCR were included to see whether lineages of non-A. hydrophila species associated strongly with MAS outbreaks (Supplementary Table 2). Five microlitres of each PCR reaction was separated on a 1.5% UltraPure Agarose-1000 gel (Invitrogen, Thermo Fisher Scientific) in chilled 0.5 × TAE buffer. A GeneRuler 100 bp Plus DNA Ladder (Thermo Fisher Scientific) was used to aid DNA profile comparisons. Following electrophoresis, the gels were stained for 30 min in 1 mg/mL ethidium bromide and de-stained in Milli-Q water (Merck Millipore, UK) for 1 h. DNA profiles resulting from the rep-PCR were visualized and gel images captured using the InGenius system (Syngene, UK). Numerical analysis of the DNA fingerprints was performed using Gel Compar II software (Applied Maths NV, Belgium). Dendrograms were constructed using the unweighted pair group method with arithmetic mean (UPGMA) and Pearson similarity coefficient. Clusters of similar banding profiles (similarity of at least 95%) were defined as an equivalent rep-PCR type where the banding patterns were considered indistinguishable by manual inspection of the respective rep-PCR gel profiles.

A subset of 12 field isolates positive by 16S rRNA PCR, representative of the main groups by rep-PCR and collected from different provinces and sampling times, were analyzed by API 20E biochemical profiling (bioMérieux, USA) and oxidase test (Oxoid) to further support the species identification (Supplementary Table 2). A. hydrophila subsp. hydrophila ATCC 7966^T was included for comparison.

A subset of 132 field isolates were selected for ngsMLST profiling such that they represented the diversity of A. hydrophila diagnostic PCR outcomes, rep-PCR profiles and sample sites. Four A. hydrophila reference isolates of varied host origin and six type strains of Aeromonas spp. were included as control material for ngsMLST library preparation (Supplementary Table 2). A total of nine MLST loci and the full-length 16S rRNA gene (Weisburg et al., 1991) were amplified by PCR, barcoded, and sequenced using the high-throughput sequencing HiSeq platform (Illumina, USA) in order to assess the genetic variation within the field isolate collection. Six loci were included (gyrB, groL, gltA, metG, ppsA, recA) based on the MLST scheme of Martino et al. (2011) with three additional markers (rpoD, dnaX, dnaJ) selected for their ability to inform phylogeny within the Aeromonas genus (Nhung et al., 2007; Martinez-Murcia et al., 2011). Primer sequences of the MLST loci and full-length 16S rRNA gene are listed in Supplementary Table 3.

The ten loci of interest were amplified individually in a 5-μL PCR containing 1 × MyTaq HS mix (Meridian Bioscience, UK), 0.2 μM of each primer pair and 2.5 ng genomic DNA template. The touchdown cycling protocol consisted of an initial denaturation at 95°C for 3 min; then 10 cycles of 95°C for 15 s, annealing at 65°C (with this temperature decreasing by 1°C each cycle), and extension at 72°C for 30 s; followed by eight PCR cycles performed at an annealing temperature of 55°C. The full-length 16S rRNA gene PCR was conducted in the same PCR master mix conditions with an initial 3 min denaturation step at 95°C; 15 cycles of 95°C for 15 s, 45°C for 30 s, and 72°C for 45 s; and a final extension at 72°C for 5 min.

Amplicons were visualized on a 1.5% TAE agarose gel and the PCR products from each isolate combined to normalize total DNA template according to band intensity and size. Samples were purified by AxyPrep Mag PCR bead clean-up kit (Axygen Biosciences, USA) at a bead ratio of 0.6 × and diluted to 0.2 ng/μL according to Qubit dsDNA HS quantification (Invitrogen; Thermo Fisher Scientific). Libraries were constructed in a miniaturized volume modified from an existing Illumina Nextera XT Library Preparation protocol (Illumina, USA), starting with 0.1 ng DNA template in a 2.5 μL tagmentation mixture and 5 μL barcoding PCR. Then libraries were cleaned up by AxyPrep magnetic beads (Axygen Biosciences, USA) at a 1:1 ratio. Finally, quantified libraries (Qubit dsDNA HS) were normalized by abundance of DNA, divided into four pools, and submitted for sequencing on a NovaSeq 6000 sequencing platform (Illumina, USA) at PE250 reads (Novogene, UK).

Mass-parallel molecular typing data (ngsMLST sequences) were filtered for quality (QC > 20), length (150 nt), and absence of primers/adaptors and complexity (entropy > 15) using fastp (Chen et al., 2018). Gene sequence assemblies and typing were performed using SRST2 v0.2.0 (Inouye et al., 2014) and PubMLST (Jolley and Maiden, 2012). The resulting sequences of nine MLST loci and the full-length 16S rRNA gene with average coverage above 400 × were concatenated and aligned using GramAlign v3.0 (Russell, 2014). A phylogenetic tree was generated by PhyML v3.3.20200621 (Guindon et al., 2010) and this included concatenated ngsMLST sequences derived from published genomes of two reference isolates (vAh ST251 strain AL09-71 and A. veronii biovar sobria LMG 13067) and 12 Aeromonas spp. type strains (Supplementary Table 2).

Fourteen isolates positive by the A. hydrophila 16S rRNA diagnostic PCR from the predominant rep-PCR types, and the reference isolate A. hydrophila T4 that originated from L. rohita in Bangladesh (Poobalane et al., 2010), were selected for WGS. Genomic DNA samples were submitted for DNA library preparation and microbial WGS (Novogene, UK) on the NovaSeq 6000 sequencing platform (Illumina, USA) at PE150 read length.

Reads from the 15 sequencing libraries were used separately during the assembly process and reads were filtered as described in Section 2.8. Raw data were assembled using Spades v3.14.0 (Bankevich et al., 2012) and the Unicycler-pipeline v0.4.8 (Wick et al., 2017). The initial de novo output was re-aligned against both A. dhakensis CIP 107500^T (assembly GCF_000820305.1), a closely related species frequently misidentified as A. hydrophila, and A. hydrophila subsp. hydrophila ATCC 7966^T (assembly GCF_000014805.1), using CONTIGuator v2.7.5 (Galardini et al., 2011) to order the contigs when continuity was incomplete. Finally, Pilon v1.23 (Walker et al., 2014) was used for polishing and correcting sequencing errors and to recover closed circular plasmid sequences. Subsequently, all genomes were annotated by DFAST v1.2.4 (Tanizawa et al., 2018) and plasmids identified using blastN v2.11.0 (Nowicki et al., 2018) against the National Center for Biotechnology Information (NCBI) plasmid database [2021–12–01]. Furthermore, genome sequences (including complete sequences, draft assemblies, and raw reads) of 116 publicly available A. dhakensis and A. hydrophila strains, and four other Aeromonas spp. to act as outgroups, were downloaded from the European Bioinformatics Institute (Supplementary Table 4). Unassembled samples were assembled according to the same process as the newly sequenced isolates.

To support species assignments, FastANI v1.33 (Jain et al., 2018) was applied to the genome sequences to calculate an ANI index against several type strains, specifically A. dhakensis CECT 7289^T and CIP 107500^T (assemblies GCF_000819705.1 and GCF_000820305.1, respectively), and A. hydrophila subsp. hydrophila ATCC 7966^T (GCF_000014805.1). An ANI value of more than 0.96 was selected as the species threshold (Ciufo et al., 2018).

MLST types were determined for each genome with mlst v2.19.0 (Seemann, 2022) and PubMLST (Jolley and Maiden, 2012). The Phyloviz v2.0a (Nascimento et al., 2017) tool was used to run the goeBURST nLV algorithm and visualize trees based on the probable patterns of evolutionary descent between allelic profiles.

PIRATE v1.0.4 (Bayliss et al., 2019) was used to create a detailed pan-genome and to classify the core and accessory genomes of all available A. dhakensis and A. hydrophila isolates. Analysis of pan-genome outputs was performed using R v4.0.0 (R Core Team, 2022). Data on clustering and presence/absence, and trees from PIRATE, were visualized using panX release bb56978 (Ding et al., 2018) and phandango v1.3.0 (Hadfield et al., 2018).

Screening for the presence of antibiotic resistance genes was conducted for the genomes with ABRicate v1.0.0 (Seemann, 2020), using multiple databases [2021–12–01]: NCBI (Feldgarden et al., 2020) and PlasmidFinder (Carattoli et al., 2014).

Core genes of the two most prevalent groups of strains observed in this present study were used to design pairs of discriminating primers, with conserved regions within core genes used to identify sites for the anchoring primers. Then, whole genomes were scanned to confirm the specificity of the primer sets using ecoPrimer v0.5 (Riaz et al., 2011). Candidate PCR primer pairs (targeting sequences in yjcS and intA_5 genes) were evaluated against a panel of field isolates and reference and type strains to confirm their specificity (Supplementary Table 2). Each PCR contained 5 ng genomic DNA in a 10 μL total volume containing MyTaq HS mix (Meridian Bioscience, UK) and 0.2 μM of each primer. Initial specificity testing was conducted using the following amplification conditions: denaturation at 95°C for 3 min; followed by 30 cycles of 95°C for 15 s, 52°C for 15 s, and 72°C for 10 s. PCR amplicons were separated through a 1% agarose gel (Meridian Bioscience, UK) containing 0.1 mg/mL ethidium bromide and visualized. A gradient PCR between 52°C and 62°C was conducted to optimize the PCR conditions and confirm the differential amplification of DNA from isolates representing the different groups of strains.

In total, 345 presumptive Aeromonas spp. colonies were recovered on isolation agar plates from P. hypophthalmus individuals exhibiting classical signs of MAS. Of these presumed A. hydrophila isolates, 58.6 (202/345) yielded the expected 103-bp amplicon in the PCR of Trakhna et al. (2009) that targets specifically the 16S rRNA gene of A. hydrophila, thus indicating the presence of this species (Supplementary Table 2). Of the nine farm sites sampled between 2018 and 2019, PCR-positive isolates suspected to be A. hydrophila were detected at seven sites, including co-isolation with PCR-negative colonies at six of these sites, indicative of the presence of other Aeromonas spp. Notably, of the 40 fish from which more than one isolate was recovered from the internal organs, in eight cases PCR-positive and negative colonies were isolated, indicating the presence of more than one Aeromonas spp. (e.g., Site 3, Fish 1; Supplementary Tables 1, 2).

A summary of the analyses performed for isolates included in this present study can be found in Supplementary Figure 1.

Analysis of the 345 field isolates by rep-PCR revealed complex DNA profiles consisting of 12–14 fragments, with each ranging from ca. 300–3,000 bp (Supplementary Figure 2). Of the 202 isolates positive by 16S rRNA PCR and thus suspected to be A. hydrophila, two distinct profiles were dominant: these were confirmed by manual inspection and numerical analysis and termed rep-PCR types A and B (Supplementary Table 2). Rep-PCR types A and B were widely distributed across the eight provinces sampled (Supplementary Figure 2), with rep-PCR type A being most prevalent (151/202; 74.8%) among the suspected A. hydrophila isolates, especially in samples collected more recently. Of note, from one fish (Site 3, Fish 3) was isolated three PCR-positive colonies, where two isolates clustered within rep-PCR type A whilst the third was rep-PCR type B (Supplementary Tables 1, 2).

More variable rep-PCR profiles were observed for the colonies isolated on Aeromonas selective medium that were negative by the 16S rRNA PCR (Supplementary Table 2, Supplementary Figure 3), and these presumed Aeromonas spp. isolates formed relatively diffuse clusters (Supplementary Figure 2) compared to the tight clusters corresponding to the two main rep-PCR types, A and B of suspected A. hydrophila isolates that tested positive by 16S rRNA PCR.

Biochemical characterisation of A. hydrophila subsp. hydrophila ATCC 7966^T and 12 field isolates, suspected to be A. hydrophila by 16S rRNA PCR and representative of rep-PCR types A and B, revealed identical API-20E biochemical profiles (704126) matching A. hydrophila, according to BacDive [https://bacdive.dsmz.de/api-test-finder].

To provide greater confidence into the identities and relatedness of the field isolates, ngsMLST was performed whereby fragments (between 477 and 1,496 bp) at each of the 10 loci were amplified for 132 representative field isolates, four A. hydrophila reference isolates from varied hosts, and six Aeromonas spp. type strains. These isolates were used as template for ngsMLST library preparation and were selected for diversity of sampling location, and the outcomes of 16S rRNA PCR and rep-PCR assays (Supplementary Table 2). After applying the selected depth threshold, a phylogenetic tree was generated for 92 libraries (87 field isolates, three reference isolates, and two type strains) based on the concatenated sequence (ca. 6,910 bp in length) from each isolate (Figure 1). Amongst the field isolates suspected to be A. hydrophila (i.e., positive by 16S rRNA PCR), two major clusters formed according to multilocus phylogenetic analysis (MLPA). Application of the PubMLST scheme and allelic profiling of the six housekeeping genes (gyrB, groL, gltA, metG, ppsA, and recA) permitted the assignment of sequence types (STs) to each of these two main clusters, namely ST656 (n = 39) and ST251 (n = 13), which corresponded to rep-PCR types A and B respectively. ST656 is associated with A. dhakensis (Jolley et al., 2010), a species closely related to A. hydrophila (Beaz-Hidalgo and Figueras, 2013), while ST251 is a hypervirulent lineage of A. hydrophila (vAh). In contrast, the PCR-negative Aeromonas spp. isolates clustered more variably and distinctly from the main A. dhakensis and vAh ST251 clades and formed a diffuse cluster that affiliated most closely to published loci of A. veronii NCIMB 13015^T (n = 26) with a minor sub-cluster related to A. jandaei CECT 4228 and A. sobria NCIMB 12065^T (n = 9; Figure 1). The remaining three field isolates from P. hypophthalmus presented as outlier groups according to MLPA. The majority of the allelic profiles of these Aeromonas spp. isolates did not match any ST in the PubMLST Aeromonas database (Supplementary Table 2).

Genome sequencing generated a mean of 12 million short-reads for each of 14 representative field isolates and the A. hydrophila T4 reference strain (Table 1). The genomes ranged between 4.82 and 4.98 Mb, with GC ratios ranging from 60.9 to 61.4%; several closed plasmids were detected consistently (Table 2). All genomes were (re-)annotated for gene location consistency using DFAST (Supplementary Table 4), aligned against selected type strains to determine ANI (Supplementary Table 5), and re-classified (Figure 2). Field isolates fell within two clades that associated with either A. dhakensis or the clonal vAh ST251 lineage within A. hydrophila, respectively. Moreover, the ANI calculations supported the existence of these two predominant species amongst the isolates collected from P. hypophthalmus (Figure 2, Supplementary Table 5). ANI values of 97.1–97.2% were estimated for nine of the genomes (all rep-PCR type A) against the two A. dhakensis type strains (CECT 7289^T and CIP 107500^T). Similarly, the five vAh ST251 genomes (all rep-PCR type B) associated most closely with the A. hydrophila subsp. hydrophila ATCC 7966^T type strain (mean ANI value of 96.8%). The reference strain, A. hydrophila T4, was also placed within the major A. hydrophila clade (ANI values of 97.0% compared to A. hydrophila subsp. hydrophila ATCC 7966^T), although separately from the ST251 lineage.

Genomes generated in this present study and all those publicly available for A. dhakensis and A. hydrophila were analyzed by MLST for the Aeromonas genus to identify STs and perform network analysis. Most of the genomes had profiles corresponding to a known ST, although 64 out of 135 analyzed genomes presented with novel gene sequences and MLST associations (denoted by “u” in Figure 3 and Supplementary Table 5). Network eBurst analysis based on sequence distance (Figure 3) revealed high diversity within and between the A. dhakensis and A. hydrophila species, thus supporting their separation into two species. For both species, strain STs were distributed typically as part of a network of single locus variants or less commonly as unique doublets or singletons (Figure 3). The A. dhakensis ST656 isolates associated with P. hypophthalmus were located within the main clonal complex of A. dhakensis strains; conversely, A. hydrophila ST251 isolated from P. hypophthalmus presented as a singleton.

Core (n = 2,517) and accessory genes (n = 33,952) were inferred to create a detailed pan-genome and to classify core and accessory genomes using all available A. dhakensis and A. hydrophila genomes (threshold of 95% presence in available genomes; Figure 4). When analyzing A. dhakensis, A. dhakensis ST656, A. hydrophila and A. hydrophila vAh ST251 pan-genomes separately, there were 2,196 shared core genes identified between the groups (threshold of 95% presence in available genomes for each group; Supplementary Figure 4). A characteristic pattern of genes within the core genomes distinguished the two closely related species of A. dhakensis (n = 1,505 unique genes not found in A. hydrophila) and A. hydrophila (n = 923 unique genes not found in A. dhakensis). Strain-specific regions associated with the A. dhakensis ST656 (n = 769 genes) and vAh ST251 (n = 425 genes) genotypes were detected when compared to all published A. dhakensis and A. hydrophila genomes, respectively (Supplementary Figure 4).

Screening for antibiotic resistance genes was conducted for each of the 135 A. dhakensis and A. hydrophila available genomes, including four Aeromonas spp. genomes as outliers (Figure 2, Supplementary Table 5). Interestingly, the genomes of the field isolates possessed sets of genes implicated in a multi-drug resistance trait and in patterns not detected in closely-related A. dhakensis and vAh ST251 counterparts (Figure 2, Supplementary Table 5). The nine field isolates of A. dhakensis ST656 each included resistance determinants to sulphonamides (sul1), trimethoprim (dfrA1), tetracycline (tetA), and a plasmid-mediated quinolone resistance gene (qnrS2). Of note, the sulphonamide (sul1) and trimethoprim (dfrA1) resistance genes were detected in both A. dhakensis ST656 and vAh ST251 outbreak genomes. Four of the vAh ST251 field isolates originating from multiple provinces in 2014 and 2015 also contained additional genetic elements conferring reduced susceptibility to gentamicin [aac(6′)-Ib4] and rifamycin (arr-2). Intriguingly, the earliest vAh ST251 isolate (TN1 collected in 2013) lacked most of these additional antibiotic-resistance determinants (Figure 2).

Gene determinants encoding reduced susceptibility to β-lactams and carbapenems were common in all A. dhakensis and A. hydrophila genomes, including in field isolates. The β-lactamase gene, bla_AQU, was found uniquely in A. dhakensis, including the ST656 outbreak strains isolated from striped catfish, while the chromosomally mediated Class D OXA β-lactamase genes, bla_OXA−726 or related bla_OXA−724, were common to genomes of both species. Sequences with homology to CphA-type carbapenemases, belonging to metallo-β-lactamases subclass B2, were also shared amongst the two species, with the cphA3 subtype associated with vAh ST251, including in the five genomes of the vAh ST251 field isolates sequenced in this present study. In contrast, a predicted cephalosporin-resistant phenotype was limited to A. hydrophila, including in the vAh ST251 field isolates (cepH) and in the A. hydrophila T4 reference genome (cepS). Markers for reduced susceptibility to chloramphenicol, aminoglycosides, macrolides and quinolones were present but distributed more sporadically across the entirety of genomes examined.

Two primer sets based on the core genes in the A. dhakensis and vAh ST251 genomes were evaluated for their ability to distinguish these two key strain groups (Supplementary Table 6). To confirm the specificity of each set of primers, the two pairs were screened in an initial evaluation against a panel of 14 Aeromonas isolates that included field isolates of A. dhakensis ST656 (n = 2) and vAh ST251 (n = 2), non-vAh A. hydrophila reference strains (n = 2), A. hydrophila subsp. hydrophila ATCC 7966^T, the A. dhakensis type strains (n = 2), and type strains of other Aeromonas spp. (n = 5). The PCR assay designed to target yjcS, which encodes a metallo-β-hydrolase in the A. dhakensis genome, yielded the predicted 223-bp amplicon at 62°C only for the DNA samples derived from the two A. dhakensis isolates collected from striped catfish in this present study (TN5 and TN49) and the two A. dhakensis type strains (CIP 107500^T and CECT 7289^T). Similarly, a PCR assay targeting the intA_5 integrase gene in the vAh ST251 genome yielded a 283-bp amplicon at the optimum 58°C annealing temperature for only the two vAh ST251 field isolates collected herein (TN1 and TN27). No amplification in either assay was observed for any non-target isolates (Supplementary Figure 5).