Setup and population

This microbiological survey study was performed at the Microbiology Laboratory of the Emek Medical Center (EMC) in Afula, Israel. The EMC laboratory serves as the regional laboratory for a population of about 0.5 million of Clalit Health Services members, living predominantly in rural settlements and small urban centers.

Whole-genome sequencing and bioinformatics analysis were performed at the Tel Aviv Sourasky Medical Center Microbiology laboratory.

Study design

The study was part of a prospective survey of bacterial gastroenteritis conducted between December 2021 until October 2022 [6]. All stool samples submitted to the EMC laboratory underwent PCR testing for Aeromonas according to the routine laboratory protocol. Samples that tested positive for Aeromonas by PCR were then cultured (details below), and one Aeromonas isolate per patient was included in the study.

The study was approved by the Ethical Review Board of the EMC.

Identification and cultivation of Aeromonas from stool samples

Stool samples were transported from the community clinics and were tested daily except on weekends, where part of each sample was suspended in ASL buffer (Qiagen, Hilden, Germany) and refrigerated at 4 °C until tested by PCR only (culturing was done from the original sample tube). Multiplex PCR for Aeromonas and other bacterial pathogens was performed as previously described [7]. Stool samples that were positive by PCR for Aeromonas were inoculated into alkaline peptone water 0.5 M NaCl with Cephalothin (10 mg/l) and incubated overnight [8] followed by sub-culturing onto SS agar plates (Hylabs, Rehovot, Israel). Presumptive Aeromonas colonies were identified by MALDI Biotyper Sirius system (Bruker Daltonics, Bremen, Germany) using the MBT IVD Library Revision software. Definite determination of species was based on whole genome sequencing as described below.

Antimicrobial susceptibility testing (AST)

AST were performed at the EMC laboratory using the Sensititre GN3F microdilution panel (Thermo Fisher Scientific) following the manufacturer’s instructions. AST included amikacin, ampicillin, ampicillin-sulbactam, aztreonam, cefazolin, cefepime, cefoxitin, ceftazidime, ceftriaxone, cefuroxime, ciprofloxacin, gentamicin, meropenem, piperacillin-tazobactam, tetracycline, tobramycin and trimethoprim-sulfamethoxazole (SXT). AST breakpoints were interpreted (when available) in accordance with the CLSI recommendations for Aeromonas species [9] or if absent (including ampicillin-sulbactam, aztreonam, cefazolin, cefepime, cefoxitin, cefuroxime and tobramycin), in accordance with the recommendations to Enterobacterales [10].

Whole genome sequencing and bioinformatic analyses

Whole genome sequencing (WGS) was done using the Illumina NextSeq500/550 system. Libraries were prepared using Illumina DNA Prep (Illumina,20,060,059). The IDT for Illumina DNA/RNA UD Indexes (Illumina, 20,027,213) were used to tagment the DNA libraries for sequencing. After sequencing of each library, FASTQ files were imported into CLC Genomics Workbench version 23.0.5(Qiagen, Denmark). Following sequencing, the reads underwent quality trimming and contigs were assembled using the CLC Genomics Workbench version 23 (Qiagen, Denmark), with the following sequences applied as a template: A. caviae NZ_AP022254, A. hydrophila CP000426, A. media CP118939, A. dhakensis CP000462 and A. veronii NZ_LKKE01000001-NZ_LKKD01000048. Identification of resistance mechanisms was conducted through a combined approach, which involved (i) annotating and identifying known acquired antibiotic-resistant genes using the Comprehensive Antibiotic Resistance Database (CARD) [11], with a cutoffs of 60% identity and 80% coverage; (ii) detection of nonsynonymous mutations of selected genes and examining their correlation with resistance traits using a multi-sequence alignment approach by AliView version 1.27. Identification of protein domain sites was done using the InterProScan software [12], in comparison with A. caviae WP_063864115. Species designation was based on the reported average nucleotide identity. This method involved comparing the genome sequence of the tested isolates with those of reference type strains accessible in the GenBank database [13].

Core genome (cg) multi locus sequence typing (MLST) scheme and phylogenetic analyses

We employed the chewBBACA [14] to develop two distinct whole-genome (wg) sequencing schemas, later refined into core genome (cg) schemas. The first schema encompassed 67 isolates representing five different Aeromonas species, and the second was exclusively tailored to the 58 isolates of A. caviae. Prodigal [15] (version 2.6.3) facilitated the identification of Coding DNA Sequences (CDS) for both schemas, while BLAST+ [16] (version 2.9.0) was used for conducting allelic comparisons. The initial phase involved constructing a wgMLST schema that included all CDS from the isolates, followed by the removal of paralogous alleles to establish the cgMLST schema. During the “CreateSchema” phase, each genome was annotated for pairwise comparisons, leading to an extensive all-against-all BLASTP search. The resulting BLAST score ratio (BSR) was calculated, with genes encoding identical or nearly identical proteins (BSR exceeding 0.6 by default) consolidated into a single database, representing alleles of the same locus.

For the A. caviae-specific schema, a cross-reference with the sequence types from PubMLST was also performed to enhance the contextual understanding of the isolates. The allelic profiles derived from the chewBBACA cgMLST schemas of both groups were then subjected to phylogenetic analysis using GrapeTree software (version 1.5.0). Neighbor Joining (NJ) and Minimum Spanning Trees were constructed for each schema, with the trees generated from the allelic profiles being visualized in iTOL [17]. This comprehensive approach allowed for a detailed exploration of the genetic relationships within and between the diverse Aeromonas spp., with a specific focus on the A. caviae isolates. Two isolates derived from two separate samples of the same patient were included in the analysis. As these two isolates were assumed to be directly related, they were included in order to provide an epidemiological control for the phylogenetic schema, i.e., the number of gene-differences that can be expected to reflect direct epidemiological linkage.

Data availability

This Whole Genome Shotgun project has been deposited in NCBI under BioProject accession number PRJNA1040111. In addition, all 58 A. caviae isolates were deposited in PubMLST under ID BIGSdb_20231207071129_2157505_94891.

Statistical analysis

MIC50 and MIC90 were calculated using RStudio software, version 4.1.2. Graphs were drawn using the ggplot2 package in R.