Taxonomic Identification of Different Species of the Genus Aeromonas by Whole-Genome Sequencing and Use of Their Species-Specific β-Lactamases as Phylogenetic Markers

Some Aeromonas species, potentially pathogenic for humans, are known to express up to three different classes of chromosomal β-lactamases, which may become hyperproduced and cause treatment failure. The aim of this study was to assess the utility of these species-specific β-lactamase genes as phylogenetic markers using whole-genome sequencing data. Core-genome alignments were generated for 36 Aeromonas genomes from seven different species and scanned for antimicrobial resistance genes. Core-genome alignment confirmed the MALDI-TOF identification of most of the isolates and re-identified an A. hydrophila isolate as A. dhakensis. Three (B, C and D) of the four Ambler classes of β-lactamase genes were found in A. sobria, A. allosacharophila, A. hydrophila and A. dhakensis (blaCphA, blaAmpC and blaOXA). A. veronii only showed class-B- and class-D-like matches (blaCphA and blaOXA), whereas those for A. media, A. rivipollensis and A. caviae were class C and D (blaCMY, blaMOX and blaOXA427). The phylogenetic tree derived from concatenated sequences of β-lactamase genes successfully clustered each species. Some isolates also had resistance to sulfonamides, quinolones and aminoglycosides. Whole-genome sequencing proved to be a useful method to identify Aeromonas at the species level, which led to the unexpected identification of A. dhakensis and A.rivipollensis and revealed the resistome of each isolate.


Introduction
The taxonomy of the genus Aeromonas is quite complex. The history of this genus of Gammaproteobacteria has been in constant change, as described by Fernandez-Bravo and Figueras [1]. The members of the Aeromonas genus are Gram-negative rods, oxidase-and catalase-positive, and glucose fermenters. They are capable of degrading nitrates to nitrites and are resistant to the vibriostatic factor O/129 (2,4-diamino-6,7-di-iso-propylpteridine phosphate) [1]. Despite these established biochemical characteristics, identifying species within the Aeromonas genus has always been a convoluted task. Included in the Vibrionaceae in 1965, the Aeromonas were subsequently classified in their own family, the Aeromonadaceae, in 1976 [2]. Since then, 36 different species have been included in the genus, up to 19 of which are considered to be potential human pathogens related to gastrointestinal, as well as blood and wound, infections in both immunocompromised and immunocompetent patients [2]. The most pervasive species are A. caviae, A. hydrophila, A. veronii and A. dhakensis [2][3][4][5]. except for two: the A. allosaccharophila reference strain (ATCC 51208), which was identified as A. veronii by MALDI-TOF, and the A. caviae D547 strain, which was identified as A. hydrophila by MALDI-TOF.
The core-genome sequence analysis allowed us to correctly identify all clinical isolates to the species level. WGS analysis also revealed which Aeromonas species share the highest core-genome similarities, thus producing three main phyllogenetic groups: one group with A. media, A. rivipollensis and A. caviae, a second one with A. dhakensis and A. hydrophila and a third one with A. sobria, A. allosaccharophila and A. veronii.
The ANI genome comparison supports the core-genome analysis, and is shown in the File S1. Additionally, type strain information is provided in Table S1.

Antimicrobial Susceptibility Testing
The Antimicrobial Susceptibility Testing (AST) was performed in twelve of the thirteen clinical isolates. A. sobria D176 was removed from the AST because it did not match the growth standards for the test. Results of AST are shown in Table 2. Maximum likelihood phylogenetic tree derived from the core-genome alignment of the Aeromonas spp. isolates and type strain genomes. All strains with names beginning with D are isolates from the laboratory. Bootstrap values above 80% are shown in branch nodes (premutation = 100). The different clusters within the phylogenetic tree contain previously identified sequences posted on GenBank. Type strain information and accession numbers are presented in Table S1. Figure 1. Maximum likelihood phylogenetic tree derived from the core-genome alignment of the Aeromonas spp. isolates and type strain genomes. All strains with names beginning with D are isolates from the laboratory. Bootstrap values above 80% are shown in branch nodes (premutation = 100). The different clusters within the phylogenetic tree contain previously identified sequences posted on GenBank. Type strain information and accession numbers are presented in Table S1. Practically all Aeromonas isolates studied were resistant to ampicillin, amoxicillinclavulanate and cefazolin, with three exceptions: A. hydrophila D173 isolate, which was only resistant to first-and second-generation cephalosporins and A. veronii (D551), susceptible to amoxicillin-clavulanate (Table 1). Only two strains showed resistance to third-generation cephalosporins and all of them were resistant to carbapenems. We also found resistance to nalidixic acid in three A. caviae and two A. veronii strains (one of them was also resistant to cotrimoxazol). No species-specific resistance patterns were observed.

β-Lactamases
Given that WGS is a technique not yet available for all clinical microbiology laboratories, and because various β-lactamases have been described in the genus Aeromonas, we hypothesized that these enzymes could help achieve a correct identification when standard laboratory techniques are inconclusive.
The assembled genomes obtained from the WGS analysis and the genomes obtained from Genbank were analysed to look for molecular mechanisms associated with β-lactam resistance using the public ResFinder and CARD databases ( Table 3). As expected, different classes of β-lactamase genes in the Aeromonas genus were found.
Thus, the phylogenetic analysis of the β-lactamase genes was carried out by UP-GMA and revealed the greatest diversity in the AmpC group (Figure 2), separating the FOX and MOX clusters. This group also possessed the highest species specificity. Following ResFinder nomenclature, all A. caviae strains were clustered in the MOX group along with the A. hydrophila strains. Nevertheless, A. dhakensis, which is phylogenetically related to A. hydrophila, showed a blaCMY-8b gene that has 84% identity with blaMOX genes. The A. dhakensis D547 and the A. hydrophila D173 strains expressed the carbapenemases CphA or Imi, which share a similarity of 94.49% and are classified in the same β-lactamase group ( Figure 2). UPGMA tree containing all β-lactamase resistance genes found in Aeromonas isolates. Class B β-lactamase branch is coloured in blue, Class C β-lactamase branch is coloured in green, and Class D β-lactamase branch is coloured in red. Analyses of the 65 nucleotide sequences were performed using MEGA X.
According to the clusters defined by WGS, the first one was composed of A. sobria, A. allosaccharophila, and A. veronii, which, despite their proximity did not share the same types of β-lactamase genes. While the blaCphA and blaAMPS genes were present in all three species, blaFOX was not detected among the A. veronii strains.
Finally, the A. media and A. caviae strains expressed the cephalosporinases blaCMY-8b or blaMOX-9 (both with a 96% identity) and the blaOXA-427 gene, which have between 96.48% and 98.36% similarity with AmpS ( Figure 2). Therefore, although each cluster was characterized by particular β-lactamase genes, the results did not reveal a sufficiently clear pattern to establish an algorithm for the identification of Aeromonas spp. It would also be important to unify a single nomenclature between the different databases consulted. UPGMA tree containing all β-lactamase resistance genes found in Aeromonas isolates. Class B β-lactamase branch is coloured in blue, Class C β-lactamase branch is coloured in green, and Class D β-lactamase branch is coloured in red. Analyses of the 65 nucleotide sequences were performed using MEGA X.
Thus, the phylogenetic analysis of the β-lactamase genes was carried out by UP-GMA and revealed the greatest diversity in the AmpC group (Figure 2), separating the FOX and MOX clusters. This group also possessed the highest species specificity. Following ResFinder nomenclature, all A. caviae strains were clustered in the MOX group along with the A. hydrophila strains. Nevertheless, A. dhakensis, which is phylogenetically related to A. hydrophila, showed a bla CMY-8b gene that has 84% identity with bla MOX genes. The A. dhakensis D547 and the A. hydrophila D173 strains expressed the carbapenemases CphA or Imi, which share a similarity of 94.49% and are classified in the same β-lactamase group ( Figure 2).
According to the clusters defined by WGS, the first one was composed of A. sobria, A. allosaccharophila, and A. veronii, which, despite their proximity did not share the same types of β-lactamase genes. While the bla CphA and bla AMPS genes were present in all three species, bla FOX was not detected among the A. veronii strains. Finally, the A. media and A. caviae strains expressed the cephalosporinases bla CMY-8b or bla MOX-9 (both with a 96% identity) and the bla OXA-427 gene, which have between 96.48% and 98.36% similarity with AmpS ( Figure 2). Therefore, although each cluster was characterized by particular β-lactamase genes, the results did not reveal a sufficiently clear pattern to establish an algorithm for the identification of Aeromonas spp. It would also be important to unify a single nomenclature between the different databases consulted.

Characterisation of Other Mechanisms of Resistance to Antimicrobials
Aside from the expected β-lactam resistances, five of the 13 isolates (two A. veronii: D551 and D553, and three A. caviae D549, D550 and D553) also showed visible resistance to nalidixic acid. The nalidixic acid resistance was associated with substitutions identified in the gyrA gene in strains D550 and D551 (Ser83Ile) and strain D553 (Ser83Arg). In strains D549 and D552, a substitution was detected in the gyrA gene (Ser83Ile) and in the parC gene (Ser80Ile). No mutations producing resistance were found in gyrB or parE, and the mutations detected were not associated with an increased quinolone resistance.
Regarding resistances to aminoglycosides, although this drug family was not included in the AST, two of the isolates expressed aminoglycoside-modifying enzymes (AME): strain D549 carried the genes aph(3 )-1b and aph(3 )-1d, which have the highest affinity for streptomycin, whereas strain D552 possessed aadA2, which hydrolyses streptomycin and spectinomycin.
Only one A. veronii strain (D551) resistant to trimethoprim-sulfamethoxazole carried the sul1 and dfrA12 genes. Query covers and identities for both genes were 100%.

Discussion
At present, identification of Aeromonas isolates to the species level remains a complex task [6,11,14]. Nevertheless, we obtained results of high accuracy using a WGS method and core-genome alignment, as have other studies [21,26,27]. This promising technique, however, is out of reach for routine use in the laboratory, and the most commonly used method, MALDI-TOF, may misidentify samples without a sufficiently comprehensive database [13,16], as occurred in our A. allosaccharophila, A. dhakensis and A. rivipollensis strains. A. rivipollensis was described for the first time in 2015, also in Catalonia, from the Ter river and, in agreement with our results, the authors describe it as a new species related with A. media [28]. To improve identification in Aeromonas, we investigated the potential of AmpC β-lactamase genes to serve as species-specific markers and, in our experience, the PCR amplification and Sanger sequencing analysis of the AmpC β-lactamase class can improve the correct identification to the species level in the genus Aeromonas.
Aeromonas spp. are known to carry several chromosomal β-lactam resistance genes belonging to Ambler classes B, C and D [18,20,[22][23][24]. By obtaining the genetic sequences of these β-lactamase genes through WGS, we hoped to shed light on whether the class C cephalosporinases MOX and FOX, now pervasive among the Enterobacteriaceae, originate from the Aeromonas genus. Ebmeyer et al. [29] detected different MOX-type enzymes in the genome of various species of Aeromonas, including A. sanaralli (MOX-1), A. caviae (MOX-2) and A. media (MOX-9). The same authors have also reported that the origin of FOX enzymes could be in the chromosome of A. allosaccharophila [30].
The high-resolution method of WGS allowed almost every Aeromonas species to be matched with a sequence type genome. The phylogenetic tree confirmed that A. sobria and A. allosaccharophila differ considerably from the genetically closest species A. veronii in the core-genome aligment, as previously described by MLPA or by gyrB and rpoB phylogenetic analysis [7,8]. The WGS approach also led to the unexpected identification of an A. dhakensis isolate, originally identified as A. hydrophila by MALDI-TOF, a species not previously reported as a clinical strain in Spain. WGS also identified an A. rivipollensis isolate previously identified as A. media.
Using WGS analysis, we were able to determine the resistance genes carried by the Aeromonas isolates. Three types of β-lactamase genes (bla CphA , bla AmpC and bla OXA ) were found in the isolates, in accordance with the species-specific patterns described in the literature. A. sobria, A. allosacharophila, A hydrophila and A. dhakensis carried all three types, A. veronii was the only clinical strain to possess only bla CphA and bla OXA , and both A. media, A. caviae and A. rivipollensis carried a bla AmpC and a bla OXA gene [18,20,[22][23][24].
It has already been reported that there is no agreement between the presence of so many genes encoding different β-lactamases (which should confer resistance to most βlactam antibiotics) with the actual resistance profile. A two-component regulator (TCR), closely related to the CreBC of Escherichia coli, has been described in Aeromonas. This TCR includes a putative mutant form of a transcription factor, the BlrA protein (related to the extended family of phosphorylation-dependent response regulators), whose gene was found immediately upstream from the blrB gene, encoding a predicted sensor kinase [25,31]. The presence of this operon prevents the expression of β-lactamase genes, and mutations in this system confer a high level of resistance to β-lactams. In some Aeromonas species (A. veronii, A. hydrophila and A. caviae), a frequency range of blrAB de-repression between 10 7 and 10 9 has been described, which increases the MICs of the β-lactams tested by 16 for ampicillin, 4 for imipenem, up to 16 for cephaloridine and up to 129 for cefotaxime [32].
Notably, the Ambler class C β-lactamases showed a pattern of high species-specificity, indicating that each Aeromonas species has a characteristic family of cephalosporinase genes (if any): bla FOX, in A. sobria, A. allosacharophila and A. hydrophila, bla AQU in A. dhakensis, bla CMY/MOX-9 in A. media and A. rivipollensis, and bla MOX in A. caviae. Therefore, plasmid-mediated AmpC genes are derived from the chromosomal AmpC genes of several members of the family Enterobacteriaceae, including Enterobacter cloacae (MIR o ACT group), Morganella morganii (DHA group), and Hafnia alvei (ACC group) [33]. As reported by Jacoby et al. [34], CMY is represented twice, as it has two quite different origins. Six current varieties (CMY-1,-8,-9,-10, -11, and -19) are related to chromosomally determined AmpC genes in Aeromonas spp., while the remainder are related to AmpC β-lactamases of Citrobacter freundii.
The chromosomal β-lactamase genes of A. caviae thus seem to be the progenitors of plasmid-mediated MOX-2 and MOX-4 and of some CMY-1-related enzymes. The plasmidmediated FOX enzymes seem to have their origin in the AmpC CAV-1 of A. media (first identified as A. punctata) [35]. Finally, A. jandaei and A. enteropelogenes also have their own chromosomal AmpC β-lactamase genes (AsbA1 and TRU-1, respectively) [4,36]. These data support that each Aeromonas species has its own chromosomal β-lactamase genes and is a potential progenitor of new emerging AmpC enzymes.
As it is mandatory to use more than one public resistance database, we found that, using the ResFinder and CARD, some entries were exclusive to either one or the other: bla AmpS and bla AmpH only appeared in ResFinder and bla AQU-1 , bla AQU-2 , bla MOX-9 and bla CepS were exclusive to CARD. These specific entries yielded matches with higher identities for our query sequences, indicating that completeness is essential for proper diagnosis and that more than one curated database for antimicrobial resistance genes should be used when screening isolates [37]. Nevertheless, regarding Aeromonas resistance genes, both ResFinder and CARD are good tools for screening, but we found that the percentage of identity for many species-specific resistance genes, such as bla AmpC of A. hydrophila and bla OXA -like of A. caviae, was too low.
Cases of fluoroquinolone resistance mediated by genetic mobile elements have been previously described in Aeromonas [4,10,27,38,39]. However, we did not find any qnr-like genes among the resistant isolates D549-D553, correlating with the findings of Ghatak et al. [27]. Additionally, we found D549-D552 to have point mutations in gyrA resulting in amino acid substitution (Ser83Ile) and parC (Ser80Ile), which are typically associated with quinoloneresistant aeromonads [40,41], as well as a mutation in gyrA: Ser83Arg in strain D553 not previously reported in Aeromonas spp.
Overall, WGS allowed the successful identification of all the Aeromonas isolates to the species level. β-lactamase genes, naturally associated with each species of the genus (a class B carbapenemase, a class C cephalosporinase and a class D penicillinase), were present in our isolates in different combinations: A. sobria, A. allosacharophila, A hydrophila and A. dhakensis carried all three types, A. veronii was the sole clinical strain to possess only a bla CphA -type and a bla OXA -type gene, whereas A. media, A. rivipollensis and A. caviae carried a bla AmpC and a bla OXA gene. To our knowledge, this is the first report in Spain of an A. dhakensis isolate, previously identified as A. hydrophila by MALDI-TOF.
Based on the results obtained, we can conclude that WGS is the most appropriate technique both to identify Aeromonas at the species level and to describe the different resistance mechanisms present. In addition, this method produces objective results that are comparable with those from any other laboratory anywhere. The only drawback for its laboratory application, at least at present, is the high cost.

Identification and Antimicrobial Susceptibility Testing
Thirteen strains were included in this study. Seven (D173-D176, D178-D180) were from the Colección Española de Cultivos Tipo (CECT) and six (547, 549-553) were isolated in our laboratory from faeces of patients with gastroenteritis. Strains were stored at −20 • C until use, and one copy at −80 • C.

DNA Extraction and Whole-Genome Sequencing
From single colonies grown in blood agar plates, the thirteen samples were transferred to an enriched growth broth (5 mL Luria Bertani) and cultured overnight at 37 • C. DNA extraction and purification was done according to the DNeasy ® UltraClean ® Microbial Kit protocol (Qiagen, Frederick, MD, USA).
The quality screening for DNA extraction was done by Invitrogen QUBIT ® 3.0 fluorometer, which determines the amount of extracted DNA (ng/µL) per sample. Concentrations that we deemed too low were concentrated using the SpeedVac at 45 • C for 10 min, thus lowering the final volume from 50 µL to 30 µL, and analyzed again. Spectrophotometry at 260 and 280 nm was also calculated to perform a quality control for DNA concentrations and purity using Nanodrop 2000/2000c(Thermofisher Scientific, Waltham, MA, USA).
Sequencing was performed by Sequentia Biotech S.L. (Barcelona, Spain) using Illumina technologies. The sequencing libraries were prepared using a Nextera XT kit (Illumina, Madrid, Spain) and sequenced on an Illumina MiSeq sequencer which generated 2 × 250 bp paired-end reads.

Whole-Genome Comparison
Complete genomes of Aeromonas were downloaded from GenBank (Table 3). Type strain genomes for ANI genome comparison and core genome analysis were downloaded from the NCBI Assembly database (available at https://www.ncbi.nlm.nih.gov/assembly/) (accessed on 9 February 2021). We used OrthoANIu standalone tool for ANI genome comparison (available at https://www.ezbiocloud.net/tools/orthoaniu) Not applicable) [49].

Antimicrobial Resistance Genes Phylogenetic Tree
Nucleotide sequences of the genes listed in Table 2 searched against the CARD and ResFinder databases were downloaded as FASTA files. We used BLAST to extract the β-lactamase genes from the genomes of our strains. MEGA X [53] was used for FASTA alignment and UPGMA tree calculation and visualization. All genes used for this analysis are indicated in Supplementary Table S2.