In our previous study, we selected 25 sampling sites in the coaster waters near Wenzhou in the East China Sea and collected surface seawater samples (Jin et al., 2021). We isolated culturable bacteria resistant to sulfonamides, tetracyclines, and quinolones from the seawater samples using antibiotic resistance plates and detected the target ARGs in these isolates using polymerase chain reaction (PCR). We isolated 1,605 bacterial strain with antibiotic resistance phenotypes, among which 51 isolates tested positive for resistance genes related to sulfonamides, tetracyclines, and quinolones, with 43 of them belonging to the genus Aeromonas (Jin et al., 2021). Additionally, we identified 39 MDR strains (Jin et al., 2021). Using PCR primers for the mcr genes, we identified 2 Aeromonas isolates containing the mcr-3 variants (Jin et al., 2021; Xu et al., 2021) (Supplementary Table S1). Preliminary 16S rDNA sequencing suggested that both isolates likely belonged to Aeromonas veronii and Aeromonas caviae. In this study, further identification of both isolates was performed using multilocus phylogenetic analysis (MLPA) based on six genes (gyrB, rpoD, dnaJ, gyrA, dnaX, and atpD) (Martinez-Murcia et al., 2011). Phylogenetic trees were constructed based on the concatenated sequences containing these six genes using the maximum likelihood method in the MEGA 7.1 software package. Subsequently, biochemical tests were conducted on the isolates using biochemical test kits (Tiangen, Beijing, China).

Antimicrobial susceptibility to 12 antibiotics belonging to seven antibiotic classes, including β-lactams [ampicillin (AMP) and ceftazidime (CAZ)], tetracyclines [oxytetracycline (OTC) and tetracycline (TET)], amphenicols [chloramphenicol (CHP) and florfenicol (FFC)], sulfonamides [sulfamethoxazole (SMZ) and trimethoprim (TMP)], quinolones [ciprofloxacin (CIP) and enrofloxacin (ENR)], aminoglycosides [gentamycin (GEN)], and polymyxin [colistin (COL)], was assessed using the disc diffusion method, with modifications derived from the methodology outlined by Huq (2020). Briefly, A. veronii 0728Q8Av and A. caviae 1029Y16Ac were cultured overnight in Mueller Hinton II (cation-adjusted)/CAMHB broth (Solarbio, Beijing, China). A 0.1 mL bacterial culture of each strain was evenly spread on a Mueller–Hinton (MH) agar plate. The sterile paper discs (7 mm in diameter) containing the antibiotics (10 μg/disc for APM, GEN, 30 μg/disc for CAZ, TET, OTC, CHP, FFC, 5 μg/disc for CIP, ENR, TMP, 300 IU/disc for COL, 250 μg/disc for SMZ) were placed on the surface of the inoculated MH agar plates. The plates were then incubated at 28°C, and the inhibition zones were measured after 18 h of incubation.

The antibiotics that induced resistance in either of the tested bacterial strains were selected to calculate the minimum inhibitory concentrations (MICs) against both strains. This was achieved using the broth microdilution method (BMD) recommended by the Clinical and Laboratory Standards Institute (CLSI) (CLSI, 2018, M100-S28). BMD is conducted using Mueller Hinton II (cation-adjusted)/CAMHB broth, a range of 2-fold dilutions of antibiotics (ranging from 0.5 to 256 μg/mL), and a bacterial inoculum density was adjusted equivalent to a 0.5 McFarland standard per well. The MIC breakpoints were interpreted according to the European Committee on Antimicrobial Susceptibility Testing (EUCAST)¹ and CLSI guidelines (CLSI, 2018, M100-S28). Escherichia coli ATCC 25922 was used as a quality control.

The genomic DNA of A. veronii 0728Q8Av and A. caviae 1029Y16Ac was extracted using cetyltrimethylammonium bromide method, following the manufacturer’s protocol BGI-PB-TQ-DNA-003A0. Subsequently, a 20 kb fragment library was constructed for each isolate, utilizing high-quality genomic DNA that met the requirements for whole-genome sequencing. Sequencing was carried out on a Pacbio Sequel II platform (BGI, Shenzhen, China). The reads were assembled using Canu (version 1.5). After multiple adjustments to achieve the optimal assembly, base correction, circularization, and plasmid library alignment were performed on the assembled sequences to obtain the final assembly results. Gene prediction was performed with Glimmer (version 3.02). The annotation of genes was accomplished through in-house pipeline on the RAST server² (Ebu et al., 2023). Acquired antimicrobial resistance genes were predicted using ResFinder 4.1³ with a threshold of percent identity set at 90% and a minimum coverage length of 60%. Genomic islands (GIs) were determined using IslandPath-DIMOB⁴ (Bertelli and Brinkman, 2018). Insertion sequences (ISs) were identified through ISfinder⁵ (Siguier et al., 2006). Transposons and integrons were predicted using TnCentral⁶ and Integron Finder,⁷ respectively (Damas et al., 2022). Comparative analysis of genome sequences was performed using Easyfig with a maximum e-value of 0.001 and identification threshold set at 98% and BRIG (with an identification threshold of 50%) (Tang et al., 2022). Virulence factors were predicted using the Virulence Factor Database (VFDB, http://www.mgc.ac.cn/VFs/) (Liu et al., 2022).

According to Jin et al. (2021), the prevalence of Aeromonas spp. with mcr genes in coastal waters is low (2 out of 51). A. veronii 0728Q8Av and A. caviae 1029Y16Ac show different levels of resistance to colistin and harbor different mcr-3 variants. To assess whether these mcr-3 variants affect resistance to colistin, we conducted functionality assays. To in silico verify mcr-3 variants in both A. veronii 0728Q8Av and A. caviae 1029Y16Ac, the sequences of mcr-3.16, mcr-3.3, and mcr-3-like, was aligned, respectively, referred to sequences of mcr-3 variants downloaded from NCBI.⁸ And a phylogenetic tree was constructed based on mcr-3 variants’ sequences using the Maximum likelihood method in the MEGA 7.1 software package.

To determine the function of mcr-3 variants, four DNA fragments, corresponding to each open reading frame (ORF) of mcr-3.16, mcr-3-like, mcr-3.16-mcr-3-like, and mcr-3.3 regions along with their 5′- and 3′-flanking regions, respectively, were amplified using the primers Q8MCR3.16-F/Q8MCR3.16-R, Q8MCR-3-like-F/Q8MCR-3-like-R, Q8MCR3.16-F/MCR-3-like-R, and Y16MCR3.3-F/Y16MCR3.3-R, and cloned into pUC19 vector by digestion of BamHI and EcoRI, respectively. The recombinant plasmids were separately transformed into E. coli DH5α. Transformants were screened on LB agar plates containing 50 μg/mL ampicillin.

Meanwhile, the ORFs of mcr-3.16, mcr-3-like, mcr-3.16-mcr-3-like, and mcr-3.3, were amplified using the primers Q8MCR3.16-BMQ-F/Q8MCR3.16-BMQ-R, Q8MCR-3-likeBMQ-F/Q8MCR-3-likeBMQ-R, Q8MCR3.16-BMQ-F/Q8MCR-3-likeBMQ-R, and Y16MCR3.3-BMQ-F/Y16MCR3.3-BMQ-R, and cloned into pET28a vector by digestion of XhoI and EcoRI, respectively. The recombinant plasmids were separately transformed into E. coli BL21. Transformants were screened using 50 μg/mL kanamycin and were induced by addition of IPTG. Colistin susceptibility of E. coli transformed with any of the recombinant plasmids was determined in triplicates.

To determine whether the mcr-3 variants in A. veronii 0728Q8Av and A. caviae 1029Y16Ac are located on the chromosome or plasmid, conjugation assays were conducted as described previously (Sun et al., 2016; Tang et al., 2024). The E. coli J53 was used as the recipient strain, and A. veronii 0728Q8Av or A. caviae 1029Y16Ac was the donor strain. The E. coli ECCNB20-2 carrying mcr-1 was used as a positive control donor (Chang et al., 2020). All bacterial strains were cultured to be in the logarithmic growth phase in LB broth at 37°C. Each of the donor bacteria, that is A. veronii 0728Q8Av, A. caviae 1029Y16Ac, and E. coli ECCNB20-2, was mixed with E. coli J53 in suitable proportion, and the bacterial mixture was transferred to filter paper on the LB agar and cultured for 8 h. Then the filter paper was suspended into a centrifuge tube containing LB broth and serially diluted. Each diluted bacterial solution (10 μL) was cultured onto screening LB agar supplemented with 2 μg/mL colistin and 100 μg/mL NaN₃, and with 100 μg/mL NaN₃, respectively. Conjugation transfer was determined as described by Tang et al. (2022).

Healthy ayu, weighing 20–25 g, were purchased from a commercial farm in Ningbo, China. Fish were raised in a recirculating system with filtered water at 20–22°C (Zhou et al., 2020). The fish were fed with pelleted dry food once a day and acclimatized to laboratory conditions for 2 weeks before experiments. Then, randomly allocate the fish into 21 tanks, each containing 10 individuals. The fish in each three tanks were intraperitoneally administered 100 μL of bacterial suspension of A. veronii 0728Q8Av and A. caviae 1029Y16Ac, at concentrations of 1 × 10⁷, 1 × 10⁸, and 1 × 10⁹ CFU/mL. Additionally, three tanks of fish were injected sterile PBS in volumes matching those of the experimental groups, serving as negative controls. All fish were kept at 20–22°C for 7 days to observe and record the morbidity and mortality rates. The degree of virulence, expressed as the 50% mean lethal dose (LD₅₀), was calculated using the Bliss method (Finney, 1985). All intervention measures are strictly carried out in accordance with the guiding principles of the Animal Experiment Ethics Committee of Ningbo University (No. 11102).

The complete genome sequences of A. veronii 0728Q8Av and A. caviae 1029Y16Ac were submitted to GenBank under Bioproject accession numbers PRJNA1070378 and PRJNA1070380, respectively.

Jin et al. (2021) identified both mcr-3.16 and mcr-3-like in strain 0728Q8Av and mcr-3.3 in strain 1029Y16Ac by DNA sequencing of PCR products. In this study, MLPA based on six house-keeping genes demonstrated that strain 0728Q8Av clustered together with selected A. veronii strains and strain 1029Y16Ac clustered together with A. caviae strains, which indicated that the two isolates belonged to A. veronii and A. caviae, respectively (Figure 1).

Regarding the biochemical characteristics, both A. veronii 0728Q8Av and A. caviae 1029Y16Ac exhibited numerous identical traits. For example, both isolates demonstrated growth at a concentration of 1–3% sodium chloride (NaCl) (w/v) but were unable to be maintained at >6% NaCl. They were both positive for acid formation when provided with glucose, sucrose, mannose, maltose, OPNG, mannitol, and saligenin, while showing negativity towards xylose, raffinose, and sorbose. Additionally, both bacterial isolates demonstrated the ability to produce oxidase, indole, and arginine dihydrolase (Table 1). However, the two bacterial isolates presented some differing characteristics. Specifically, A. veronii 0728Q8Av, unlike A. caviae 1029Y16Ac, tested positive in the Voges-Proskauer test and exhibited urea hydrolysis. A. caviae 1029Y16Ac could utilize arabinose and lactose, producing acid, whereas A. veronii 0728Q8Av could not (Table 1).

The susceptibility testing showed that A. veronii 0728Q8Av exhibited resistance to COL (with a MIC of 8 μg/mL), OTC (64 μg/mL), TET (32 μg/mL), TMP (32 μg/mL), SMZ (256 μg/mL), AMP (64 μg/mL), CHP (16 μg/mL), and FFC (64 μg/mL), and intermediate to GEN (8 μg/mL), CIP (2 μg/mL), and ENR (2 μg/mL), but remained susceptible to CAZ (4 μg/mL) (Table 2). In contrast, A. caviae 1029Y16Ac exhibited resistance to COL (with a MIC of 4 μg/mL), AMP (32 μg/mL), TET (32 μg/mL), OTC (32 μg/mL), TMP (256 μg/mL), SMZ (256 μg/mL), CIP (4 μg/mL), ENR (4 μg/mL), and intermediate to CHP (8 μg/mL), FFC (16 μg/mL), and GEN (8 μg/mL), but remained susceptible to CAZ (8 μg/mL) (Table 2).

A. veronii 0728Q8Av harbored one circular chromosomal genome of 4,752,946 bp with G + C content of 58.5% and no plasmid was found. A. caviae 1029Y16Ac harbored one circular chromosomal genome of 4,779,291 bp with G + C content of 60.86% and two plasmids, 1029Y16AcP1 (21,050 bp) and 1029Y16AcP2 (8,340 bp).

Acquired ARGs were identified by Resfinder 4.1 analysis. The chromosomal genome of A. veronii 0728Q8Av contained eight types of acquired ARGs indicating resistance to aminoglycoside [aac(6′)-Ib3, aph(3″)-Ib, aph(6)-Id, aadA1], polymyxin (mcr-3.16, mcr-3-like), quinolone (qnrVC4), sulfonamides (sul1, dfrA14), beta-lactam (bla_OXA-10, ampS, cphA4), tetracycline [tet(E)], amphenicol (floR, cmlA1), and quaternary ammonium compound (qacE∆1). The genome contained 15 GIs and almost all the ARGs were located on the GI1 of A. veronii 0728Q8Av (AvGI1), except for mcr-3.16, mcr-3-like, ampS and cphA4 (Supplementary Figure S1A; Supplementary Table S3). The chromosomal genome of A. caviae 1029Y16Ac contained seven types of acquired ARGs indicating resistance to aminoglycoside [aac(6′)-Ib3, aadA1], polymyxin (mcr-3.3), sulfonamides (sul1), beta-lactam (bla_VEB-3, bla_MOX-6, bla_OXA-10), tetracycline [tet(E)], amphenicol (catB3) and quaternary ammonium compound (qacE∆1). No ARGs were identified on plasmid 1029Y16AcP1 and 1029Y16AcP2. The chromosomal genome presented 23 GIs, and almost all the ARGs were located on the GI11 of A. caviae 1029Y16Ac (AcGI11), except for mcr-3.3, bla_MOX-6, and tet(E) (Table 2; Supplementary Figure S1B; Supplementary Table S3).

The virulence factors’ prediction suggested that virulence factors were significantly different between A. veronii 0728Q8Av and A. caviae 1029Y16Ac (Table 3). The Tap type IV pili of both strains exhibited similar structural features. However, A. veronii 0728Q8Av possessed 16 genes related to MSHA type IV pili, whereas A. caviae 1029Y16Ac had only 3. And A. veronii 0728Q8Av had type I and Flp type IV pili, whereas A. caviae1029Y16Ac did not. Flagella analysis revealed that both strains had polar flagella. However, A. veronii 0728Q8Av possessed 54 genes related to polar flagella, whereas A. caviae 1029Y16Ac had only 7. Both stains had a type II secretion system (T2SS) containing 14 genes. Toxin analysis revealed that both strains contained Hemolysin and Thermostable hemolysin (TH). Interestingly, A. veronii 0728Q8Av had Aerolysin, but A. caviae 1029Y16Ac did not.

Blast analysis revealed that both mcr-3.16 (1623-bp of ORF) of A. veronii 0728Q8Av and mcr-3.3 (1623-bp of ORF) of A. caviae 1029Y16Ac exhibited 100% identity to previous reported mcr-3.16 (WP-111273845.1) and mcr-3.3 (WP-099982814.1). MCR-3-like of A. veronii 0728Q8Av provided a 145-aa substitution by 29 amino acids at the carbon terminus, while a single amino acid difference occurred at the 273rd residue within the 396 amino acid residues at the Nitrogen terminus, compared to reported MCR-3-like (MF495680) (Figure 2). The phylogenetic tree unveiled a close clustering of both MCR-3.16 and MCR-3.3 with other predetermined MCR-3 variants. Additionally, the MCR-3-like variant in this study exhibited clustering alongside the one reported in A. veronii isolated from retail chicken (MF495680) (Supplementary Figure S2).

To determine the function in mediating colistin resistance of mcr-3 variants, 2,685 bp, 1,920 bp, 4,008 bp, and 2,787 bp DNA fragments, corresponding to the ORFs of mcr-3.16, mcr-3-like, mcr-3.16-mcr-3-like, and mcr-3.3 along with their 5′- and 3′-flanking regions, respectively, were amplified to construct the recombinants pUC19-mcr-3.16, pUC19-mcr-3-like, pUC19-mcr-3.16-mcr-3-like, and pUC19-mcr-3.3, and transformed into E. coli DH5α. The results showed that transformants containing pUC19-mcr-3.16, pUC19-mcr-3.16-mcr-3-like, and pUC19-mcr-3.3, had colistin MICs of 4 μg/mL, 4 μg/mL, and 2 μg/mL, respectively, which were 8-, 8- and 4-fold higher than the MIC of DH5α containing pUC19 alone (0.5 μg/mL). Transformants containing pUC19-mcr-3-like had a colistin MIC of 0.5 μg/mL, similar to DH5α containing pUC19 alone (Table 4). Meanwhile, the ORFs of the three mcr-3 variants were cloned to construct the recombinants pET28a-mcr-3.16, pET28a-mcr-3-like, pET28a-mcr-3.16-mcr-3-like, and pET28a-mcr-3.3. The results showed that the transformants containing pET28a-mcr-3.16, pET28a-mcr-3.16-mcr-3-like, and pET28a-mcr-3.3, had colistin MICs of 16 μg/mL, which were 2-fold higher than the MIC of BL21 containing pET28a alone (8 μg/mL). Transformants containing pET28a-mcr-3-like had a colistin MIC of 8 μg/mL, similar to BL21 containing pET28a alone (Table 4).

A segment consisting of mcr-3.16-mcr-3-like-dgkA was observed in A. veronii 0728Q8Av, which was identical to that in A. salmonicida Z5-5 (GCA_003265515.1), a foodborne isolate from chicken meat in China. A similar genetic composition pattern mcr-3.3-mcr-3-like-dgkA was identified in many other A. veronii isolates (MF495680, CP040717, GCA_003265545.1, and GCA_003265585.1). The mcr-3-like gene in A. veronii 0728Q8Av was located 66 bp downstream of the mcr-3.16 gene, and the dgkA gene was located 118 bp downstream of mcr-3-like, which were also found in A. veronii 172 (MF495680), A. veronii Z2-7 (GCA_003265545.1) and A. veronii ZJ12-3 (GCA_003265585.1). In comparison, a 66 bp intergenic region sequence existed between mcr-3.16 (mcr-3.3) and mcr-3-like in A. salmonicida Z5-5 (GCA_003265515.1) and A. veronii HX3 (CP040717), while the intergenic region sequence between mcr-3-like and dgkA spaned 206 nucleotides. Similar to A. salmonicida Z5-5, the 5′ flanking region of the segment mcr-3.16-mcr-3-like-dgkA in A. veronii 0728Q8Av contained an insertion sequence ISKpn3. However, in A. salmonicida Z5-5, two hypothetical proteins were inserted between mcr-3.16 and ISKpn3, distinguishing it from A. veronii 0728Q8Av. Additionally, an insertion sequence ISAs29 could be found adjacent to the dgkA gene in the 3′ flanking region of mcr-3.16-mcr-3-like-dgkA (Figure 3).

A segment consisting of mcr-3.3-mcr-3-like was discovered in A. caviae 1029Y16Ac. But the mcr-3-like gene was incomplete and only had 116 amino acids, which showed 100% identity to the Nitrogen terminus of the MCR-3-like (MF495680). Compared to other Aeromonas spp., it seemed that the 3′ end of the mcr-3-like and its downstream dgkA gene were completely absent in A. caviae 1029Y16Ac. None ORFs but 10 ISs, namely ISAs12, ISAs15, ISAs2, ISAs29, IS5D, IS50R, and ISAeme2, were observed at the flanking regions of the segment mcr-3.3-mcr-3-like in A. caviae 1029Y16Ac (Figure 3). No transconjugants were observed on the colistin-containing screening LB agar plates coated with bacterial culture from the conjugation assay using A. veronii 0728Q8Av or A. caviae 1029Y16Ac. This indicated that the mcr-3 variants, that is mcr-3.16 and mcr-3-like in A. veronii 0728Q8Av and mcr-3.3 in A. caviae 1029Y16Ac, were not transferred into E. coli J53 through conjugation, suggesting that these mcr-3 variants were located on the chromosomes of both A. veronii 0728Q8Av and A. caviae 1029Y16Ac (Supplementary Figure S3).

The results showed that gene island AvGI1 in A. veronii 0728Q8Av harbored 12 ARGs, including tet(E), aph (3″)-Ib, aph(6)-Id, floR, qnrVC4, cmlA1, bla_OXA-10, aac(6′)-Ib3, aadA1, dfrA14, truncated qacE (qacE∆1), and sul, which potentially mediated the resistance against tetracyclines, amphenicol, aminoglycosides, quinolones, beta lactams, quaternary ammonium compounds, and sulfonamides antibiotics (Figure 4A). Meanwhile, AvGI1 contained five complete ISs (ISAs1, ISVsa3, IS26, IS6100, and IS3000), one incomplete IS (IS15DI), and one transposon Tn5393. In addition, AvGI1 had a typical class 1 integron containing the integrase intI1 localized at the attI1 site and the genes qacE∆1 and sul1 conferring resistance to quaternary ammonium compounds and sulfonamides, respectively. The gene cassette array featured four attC sites, partitioning it into four gene boxes: qnrVC4, orf-clmlA, bla_OXA-10-aac (6′)-Ib3, and orf-aadA1 (Figure 4A).

The gene island AcGI11 in A. caviae 1029Y16Ac harbored 7 ARGs, including bla_VEB-3, sul1, qacE∆1, aac (6′)-Ib3, aadA1, bla_OXA-10 and catB3, which potentially mediated the resistance against beta lactams, sulfonamides, quaternary ammonium compounds, amphenicol, and aminoglycosides (Figure 4B). Meanwhile, AcGI11 contained six complete ISs, including IS6100, ISAeme2, ISAeme4, two copies of ISKpn9, and IS91, and one transposon TnAs1. Similarly, AcGI11 had the typical class 1 integron and the gene cassette array featured four attC sites, partitioning it into four gene boxes: catB3, bla_OXA-10-orf-aadA1, aac(6′)-Ib3, and catB3 (Figure 4B).

Blast analysis indicated that three genomic segments were similar to AvGI1 (query coverage >72%, identities >90%), that is A. hydrophila NUITM-VA1 (AP025277), A. caviae WCW1-2 (CP039832) and A. simaiae A6 (CP040449), while no genomic segment was similar to AcGI11. Multiple sequence alignment revealed a high similarity in the class 1 integron, including the ORFs and ISs within the gene cassettes, between A. veronii 0728Q8Av (PRJNA1070380) and A. hydrophila NUITM-VA1 (AP025277), as well as A. caviae WCW1-2 (CP039832). The transposon Tn5393 and the insertion sequence IS3000 were observed upstream and downstream of the class 1 integrons in A. veronii 0728Q8Av, A. caviae WCW1-2 (CP039832), and A. simaiae A6 (CP040449). But, A. veronii 0728Q8Av contained an additional segment, IS15DI-orf-floR-orf-ISVsa3-IS26, located between Tn5393 and the class 1 integron. A. hydrophila NUITM-VA1 (AP025277) did not form the transposon Tn5393 (Figure 5).

The mortality of fish infected with the A. veronii 0728Q8Av and A. caviae 1029Y16Ac isolate are summarized in Table 5. Using the calculated method reported by Mittal et al. (1980), the LD₅₀ value in a 7-day period obtained from A. veronii 0728Q8Av was measured at 2.15 × 10⁷ CFU/mL and was categorized as virulent, while the LD₅₀ of A. caviae 1029Y16Ac was calculated at 1.35 × 10⁹ CFU/mL and was non-virulent (Table 5).