Genomic analysis of Elizabethkingia species from aquatic environments: Evidence for potential clinical transmission

Highlights • Identification of closely related (< 50 SNV) clinical and environmental aquatic Elizabethkingia anophelis isolates.• Identification of a provisional novel species Elizabethkingia umaracha.• Novel blaGOB and blaB carbapenemases and extended spectrum β-lactamase blaCME alleles identified in Elizabethkingia spp.• Analysis of the global phylogeny and pangenome of Elizabethkingia spp.• Identification of novel ICE elements carrying uncharacterised genetic cargo in 67 / 94 (71.3%) of the aquatic environments Elizabethkingia spp.


Introduction
The environment is a known reservoir for both opportunistic pathogens and antimicrobial resistant (AMR) bacteria (Ishii, 2019). It is important to investigate environmental microbial populations as prominent extended-spectrum β-lactamases (ESBLs), quinolone resistance genes and carbapenemases originated from marine and soil bacterium and have subsequently entered clinical isolates through horizontal gene transfer (HGT) (Poirel et al., 2005;Wyres and Holt, 2018), as well as identifying potential infection transmission pathways (Ishii, 2019).
Interest in Elizabethkingia spp. is rising as they constitute difficult to treat emerging pathogens within hospitals and healthcare settings (Green et al., 2008;Lau et al., 2016;Teo et al., 2013). Elizabethkingia infections occur most often in newborns and immunocompromised patients, and the most common presentation is septicaemia (Sarma et al., 2011). However, reports of meningitis (caused by E. meningoseptica and E. anophelis), pneumonia, urinary tract infection, skin and soft tissue infections are also common (Lin et al., 2009;Venkatesh et al., 2018). Case fatality rates for Elizabethkingia spp. are high at ~25.2% in all species (Seong et al., 2020) and higher in cases of septicaemia and meningitis at 54% for E. meningoseptica (Moore et al., 2016) and 28.4% for E. anophelis infections (Lin et al., 2018b). Elizabethkingia pathogenesis is largely unknown, though several virulence factor homologs have been reported, including capsule proteins, adhesins, iron uptake proteins and proteins contributing to biofilm formation Janda and Lopez, 2017;Li et al., 2015;Liang et al., 2019).
Treatment of Elizabethkingia infections is complicated because most species are intrinsically resistant to clinically important antibiotics, including carbapenems and other β-lactams and aminoglycosides (Janda and Lopez, 2017;Li et al., 2015;Lin et al., 2019a). Resistance to fluoroquinolones and sulfamethoxazole-trimethoprim have also been observed (Lin et al., 2018a). Consistently, Elizabethkingia species genomes sequenced to date harbor multiple chromosomal antimicrobial resistance genes (ARGs), including genes bla B and bla GOB , associated with resistance to carbapenems, and extended-spectrum β-lactamase gene bla CME , conferring resistance to all cephalosporins (González and Vila, 2012). For fluoroquinolones resistance, mutations within conserved regions of DNA gyrase subunit A (GyrA) have been observed, including Ser83Ile and Ser83Arg (Jian et al., 2018). ARGs have also been identified in Elizabethkingia integrative and conjugative elements (ICEs) (Xu et al., 2019), mobile genetic elements (MGEs) capable of integrating into a host genome and propagated during chromosomal replication and cell division (Wozniak and Waldor, 2010). Only two plasmids have been described and sequenced from two Elizabethkingia species: E. anophelis strain F3201 (Xu et al., 2019) and E. miricola strain EM_CHUV (Opota et al., 2017).
Cases of Elizabethkingia infections have been increasing over the past few decades, with reports surfacing in 25 countries across six continents (Breurec et al., 2016;Burnard et al., 2020;Hsu et al., 2011;Lau et al., 2015;Lin et al., 2019b;Perrin et al., 2017;Reed et al., 2020;Teo et al., 2013). Elizabethkingia mode of transmission remains unclear, but exposure to contaminated environments, especially waterbodies, medical devices, hemodialysis and mechanical ventilation equipment, hospital fomites, water faucets and healthcare worker hands, have all been implicated (González and Vila, 2012;Lin et al., 2018a). In addition, infections caused by E. anophelis have been linked to transmission events associated with mosquitoes in the Central African Republic (Frank et al., 2013); however, this hypothesis is contentious due to a report of vertical transmission from mother to infant (Lau et al., 2015).
Elizabethkingia species are multidrug-resistant emerging pathogens with high case-fatality rates. To date, most studies on Elizabethkingia have characterized clinical isolates (Eriksen et al., 2017). However, given that waterbodies are reservoirs and implicated in Elizabethkingia transmission pathways, this study provides a genomic analysis of whole-genome sequences (WGS) derived from 94 Elizabethkingia isolates originating from aquatic environments in South Australia, as well as characterizing and comparing their genomic and antimicrobial resistance profiles to other publicly available Elizabethkingia environmental and clinical strains.

Sample collection and bacterial isolation
Water samples (~10 L) were collected monthly in triplicate from July 2018 to July 2019 in South Australia. The four locations represented two sources: (i) stagnant water (site A) and (ii) inland wetlands recharged by seasonal rainfall/runoff inflows (site B and C) or by a river (site D). Site A was a small rural reservoir created by damming a natural rainwater catchment area. It was fenced and not accessible to or impacted by livestock but was regularly visited by birds and particularly by ducks. Site B (wetland) was a recreation reserve covering an area of 19.4 hectares. Site C (wetland) covered 172 hectares and had a maximum capacity of 1200 megaliters. The distance between the two wetland sites was ~10 km. Site D was a river 2508 km in length, and the area sampled was flowing water near a wetland that covers a total of 42 hectares. Site D was used for recreational purpose only. At all sites, surface water samples were collected by dipping three sterile 10 L collection tanks below the surface. All samples were stored on ice directly after collection and processed within 2-3 h.

Isolation of carbapenem-resistant Elizabethkingia spp
Samples were processed on the day of collection. First, water samples were serially diluted and 500 μl from 2 to 3 consecutive 10-fold serial dilutions were plated in triplicate on Oxoid Brilliance™ CRE Agar plates (Thermo Fisher Scientific Australia, Adelaide, SA). Cultures were incubated at 25 • C, 37 • C and 44 • C for 24 h. Next, using pre-sterilized toothpicks, single colonies growing on CRE Agar were plated on Plate Counting Agar (PCA; Thermo Fisher Scientific). PCA cultures were incubated at 37 • C for 18-24 h, or until sufficient bacterial growth had occurred. A total of 667 bacteria were isolated and identified with Matrix-Assisted Laser Desorption Ionization-Time of Flight Mass Spectrometry (MALDI-TOF MS) (Bruker Daltonics) and preserved in glycerol stocks (40% v/v) at − 80 • C.

MALDI-TOF MS species identification
Fresh bacterial single isolates (<24 h old) were resuspended in 1 ml 70% ethanol, vortexed for 1 min, and centrifuged at 13,000 rpm for 2 min. The supernatant was removed, and the pellet re-dissolved with 5 µl of 70% formic acid (Baker; 90% stock) and 5 µl acetonitrile (CAN, LC-MS Grade, Merck). After 2 min of centrifugation at 13,000 rpm, 1 µl of supernatant was spotted onto the target plate and left to dry. The sample was overlaid with α-cyano-4-hydroxycinnamic acid (HCCA) (1 μl) matrix (10 mg/ml − 1 ) and allowed to crystallize at room temperature. One μl of Bacterial test standard (Bruker Daltonics) in 50% (v/v) ACN containing 2.5% (v/v) trifluoroacetic acid (LC-MS Grade; Thermo Fisher Scientific) was spotted, left to dry and overlaid with HCCA for calibration. MALDI-TOF MS analysis was acquired on an autoflexTM speed MALDI-TOF/TOF mass spectrometer (Bruker Daltonics) operated in linear positive mode under MALDI Biotyper 3.0 Real-time Classification (version 3.1, Bruker Daltonics) and FlexControl (version 3.4, Bruker Daltonics) software. Spectra were acquired in the mass range of 2000 to 20,000 Da with variable laser power, and a total of 1200 sum spectra were collected in 40 shot steps. The sample spectra were identified against an MSP database library (5989 MSP entries). Identification scores of 2.300-3.000 indicate highly probable species identification, scores of 2.000-2.299 indicate secure genus identification and probable species identification, scores of 1.700-1.999 indicate probable genus identification, and a score of ≤ 1.699 indicates that the identification is not reliable.

DNA extraction
Water samples were concentrated by vacuum filtration through a 0.2 μm nitrocellulose filter (Millipore) then stored at − 80 • C until DNA extraction. Total genomic DNA from each 0.2 μm filter was extracted using the DNeasy PowerWater kit (Qiagen) according to the manufacturer's instructions. DNA from MALDI-TOF MS identified colonies with scores 2.000-3.000 were extracted using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer's instructions. Nucleic acid quality (i.e., 260/280 ratio) was measured with Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific). DNA concentrations for all samples were measured by fluorometric quantitation using a Qubit instrument and High Sensitivity dsDNA HS Assay kit (Thermo Fisher Scientific), and purified DNA extracts were stored at − 20 • C until used.

Absolute quantification of Elizabethkingia spp., E. anophelis and E. meningoseptica
Standard curves to determine the absolute quantity, efficiency, linear range, and reproducibility of Elizabethkingia spp., E. anophelis and E. meningoseptica assays were prepared using the American Type Culture Collection (ATCC) strain E. meningoseptica ATCC 13,253 and the clinical isolate identified by MALDI-TOF and whole genome sequencing E. anophelis DSM 23,781. The ATCC strain and the clinical isolate were purified using the QIAamp DNA Mini and Blood Mini kit according to the manufacturer's protocol for isolation of genomic DNA from bacterial plate cultures (page 56, Handbook 05/2016, Qiagen, Sydney, NSW). The DNA concentration and quality were determined using a Qubit 2.0 fluorometer (Life Technologies). Standards were prepared by serial diluting the DNAs and by calculating E. anophelis and E. meningoseptica copy number with the following equation: Elizabethkingia copy number = (concentration of template in ng × NL) / (n × 10 9 × 660) where NL is the Avogadro constant (6.02 × 10 23 ), n is the genome length of the standard in base pairs or nucleotides and 660 is the average molecular weight of double-stranded DNA.
Digital droplet PCR (ddPCR) was used to quantify the copy numbers of the standard's serial dilutions. ddPCR was performed using QX200™ ddPCR™ Supermix for Probes (No dUTP, Biorad, Australia) and a QX200™ Droplet Digital™ PCR System with automated droplet generation (Bio-Rad, Pleasanton, CA, USA). All ddPCR amplifications were conducted in 20 μL reaction mixtures containing 10 μL of Probe Supermix, 1 μL of each individual primer (100 nM), 2 μL of template DNA, and 6 μL of ultrapure PCR-grade water. The ddPCR amplification conditions were as follows: 25 • C for 3 min, 95 • C for 10 min, 40 cycles at 94 • C for 30 s and 60 • C for 1 min, 98 • C for 10 min and 4 • C for hold.
All qPCR analyses were carried out in duplicate on a LightCycler® 480 Instrument II (Roche Life Science) with positive, negative, and nontemplate controls included. Individual real-time qPCR assays were used to quantify E. anophelis, E. meningoseptica and Elizabethkingia spp. genome copies using a multiplex probe assay with the primers and probes described in Table 1.
Amplification was done in 25 µl reaction volumes consisting of 10 µl of the LightCycler® 480 SYBR Green I Master (Roche Life Science), 5 µl of DNAse/RNAse free water (Roche Life Science), 5 µl of the primerprobe mixture, and 5 µl of template DNA within the concentration range of 40 to 50 ng/µl. The cycling conditions were as follows: initial denaturation at 95 • C for 3 min, followed by 40 cycles of 95 • C for 10 s and 64 • C for 30 s (Kelly et al., 2019). Fluorescence data were acquired at the end of the annealing step of each cycle. All mixes were made using a Biomek Automated Liquid Handler (Beckman Coulter) to avoid pipetting errors. The efficiency of the different real-time PCRs ranged from 97 to 100%. Secondary structures were not encountered in any of the runs. The threshold of each single run was placed above any baseline activity and within the exponential increase phase. The cycle thresholds (C T ) were determined by a mathematical analysis of the resulting curve using the software manufactured by Roche Life Science. The C T values of the non-template controls were always 40 or above, indicating no amplification. Dissociation curves were determined for qPCR products to confirm product integrity and the absence of PCR inhibitors. Among the different qPCR coefficients, attention was given to the R coefficient, which was used to analyze the standard curves obtained by linear regression analysis. Most of the samples, and all standards, were assessed with a minimum of two runs to confirm the reproducibility of the quantification.
Real-Time PCR datasets were analyzed using analysis of variance (ANOVA). To evaluate the absolute abundance of gene copy numbers in water samples, F-tests were used to compare variance. Normality was tested with a Shapiro-Wilks test and by inspection of residuals, and variance homogeneity by Levene's test. When data failed to satisfy one of these tests, an appropriate transformation was applied (log or squareroot transformation). Tukey's honestly significant difference (HSD) method and the modified version for unequal sample size (Unequal N HSD) were used for post hoc comparisons with a 0.05 grouping baseline. Graphs were drawn using GraphPad Prism version 9 (GraphPad Software Inc.).

Whole-genome sequencing
Whole-genome sequencing was performed as described previously (Foster-Nyarko et al., 2020). Briefly, WGS was performed on the Illumina NextSeq 500 platform using a modified Nextera low input tagmentation approach. Genomic DNA was normalized to 0.5 ng µl− 1 with 10 mM Tris-HCl before the library preparation. The pooled library was run at a final concentration of 1.8 pM on a mid-output flow cell following Illumina recommended denaturation and loading parameters. Data was uploaded to Basespace (www.basespace.illumina.com), where the raw data was converted to FASTQ files for each sample.

Minimum inhibitory concentration (MIC) testing
Representative isolates from each Elizabethkingia clade and isolates harbouring unique combinations of bla B and bla GOB genes were selected for MIC testing (n = 10) against 38 clinically relevant antimicrobials as described previously (Burnard et al., 2020). Antibiotic testing plates were hand prepared, inoculated and incubated in accordance to AS ISO 20,776.1-2017. Quality Control of antibiotic and testing isolates was in accordance to Clinical and Laboratory Standards Institute (CLSI) M100 ED31:2021; plate reading in accordance to the European Committee on Antimicrobial Susceptibility Testing (EUCAST) reading guide v 3.0 2021.
Both the guidelines of the EUCAST pharmacokinetic-pharmacodynamic (PK-PD) "non-species" breakpoints (Kahlmeter et al., 2006) and the non-Enterobacteriaceae breakpoints of the CLSI (Jorgensen et al., 2007) were used in the AMR phenotypic analysis.

Results
In this study, 94 Elizabethkingia isolates were collected from aquatic environments in South Australia from 2018 to 2019. Strains sourced from wetlands (site B & C) constituted the majority [n = 70 (B = 50, C = 20); 75%], followed by dam (site A, n = 22; 23%) and then river (site D, n = 2; 2%) samples. Associated metadata on all isolates used, including 54 sourced from outside this collection used in phylogenetic and gene screening analyses, is available in Supplementary Data 1.

Genome assembly
Draft genomes were assembled using shovill v1.0.4. Genome size ranged from 4,039,979 bp to 4,660,922 bp, with an average size of 4,459,168 bp. The number of contigs per genome ranged from 25 to 160, with a mean of 55. Read depth ranged from 23.26 to 80.63, with a mean of 38.79. Full assembly statistics can be viewed in Supplementary Data 2.

Absolute quantification of Elizabethkingia spp., E. anophelis and E. meningoseptica
Quantitative data targeting a generic Elizabethkingia spp. gene marker and E. anophelis and E. meningoseptica markers were used to estimate the absolute abundance of each in samples from the four aquatic sites. The absolute abundance of Elizabethkingia spp. in the dam sample was on average 7.6 × 10 3 gene copies/mL, representing 1.36×10 − 6 of the total bacterial community (16S rRNA qPCR-based). In the wetland samples, Elizabethkingia spp. ranged from 3.5 × 10 4 genes/mL to 4.6 × 10 4 genes/mL, representing 6.25×10 − 6 to 8.21×10 − 6 of the total bacterial community (Fig. 1). In each case, the absolute abundance of E. anophelis was a factor of ten lower than the total Elizabethkingia spp. absolute abundance, indicating that it is not the dominant species within the aquatic environments. E. meningoseptica were detected and gene copies/mL were ranging from 23 (site A; dam) to 50 copies/mL (site B, C and D; wetlands) on average.

Phylogenetic analysis
A phylogenetic tree comprised of 148 Elizabethkingia isolates was constructed using Phylosift ( Fig. 2) with 94 isolates from the Australian aquatic environments (this collection), 27 isolates from Australian clinical samples and Australian hospital environments, and 27 international strains available from Genbank. Where metadata was available, Elizabethkingia isolates were derived from the environment (n = 102), humans (n = 42), Anopheles gambiae (n = 2), and one isolate each from Zea mays (corn) and a frog. The distribution of species in the phylogenetic tree were E. anophelis (n = 52), E. meningoseptica (n = 5), E. miricola (n = 78), E. bruuniana (n = 3), E. ursingii (n = 2), E. occulta (n = 1) and a novel clade of Elizabethkingia spp. (n = 7). The seven species of Elizabethkingia were clearly separated from each other, with E. meningoseptica appearing the most distant from the other species.
The Elizabethkingia isolates from our aquatic environment study formed five clades that branched alongside international and clinical isolates. Clade 1 contained three E. anophelis from the aquatic environmental study (ER-QUAD-EK_54, QUAD-EK_14 and QUAD-EK_22) that were closely related to Australian clinical isolates (EkS2, EkQ5 and EkQ17) with an average of 36 SNV (Single nucleotide variants) across 87% of the core genome (EkS2 as reference) and hospital environment isolates EK2 and EK6 with an average of 42 SNV. The two E. anophelis isolates from dam samples, QUAD-EK_14 and QUAD-EK_22, were separated by 2 SNV, while ER-QUAD-EK_54 from the Australian wetland samples differed by an average of 33 SNV from dam isolates. Clade 2 consists of 13 clonal E. anophelis from our study that differ by an average of 2 SNV between each other and 807 SNV from our E. anophelis isolates situated in Clade 1 (pairwise SNP matrices for E. anophelis isolates provided in Supplementary Data 4). Clade 3 appears as a novel clade, represented by seven isolates (ER-QUAD-EK_21, QUAD-EK_08, QUAD-EK_09, QUAD-EK_10, QUAD-EK_16, QUAD-EK_07, and QUAD-EK_05) isolated from an Australian dam. Isolates within this apparently novel clade of the Elizabethkingia were most closely related to E. bruuniana but appeared genetically distinct in a progressiveMauve analysis (Supplementary Data 3) and differed by an average of 124,216 SNPs to E. bruuniana isolate EkQ11. Clade 4 constitutes ten E. miricola isolates with an average of 66 SNV between each other (range 0 -197 SNPs) across 83% of the core genome (EkQ1 as reference). These isolates branch alongside three E. miricola isolates from Australian clinical samples (EkQ10, EkQ13 and EkQ1) however, the SNV between these two branches is ~21,539. Clade 5 represents a group of 61 clonal E. miricola isolates (average 7 SNV) from Australian wetlands, with the closest relative strain CP03929, from a water sample from Taiwan, at ~21,629 SNPs difference. Pairwise SNP matrices for E. miricola isolates provided in Supplementary Data 4.

Identification of proposed new species Elizabethkingia umeracha sp. nov
The seven isolates in clade 3, with an average 124,216 SNV to E. bruuniana isolate EkQ11, were investigated to determine whether they constituted a closely related, yet distinct species to E. bruuniana. For this purpose, 16S rRNA and rpoB sequence identities as well as ANIb, ANIm and GGDC values (the latter mimicking DDH values) were calculated (averages presented in Table 3; full analysis in Supplementary Data 5). Except for a single ER-QUAD-EK_05 16S rRNA sequence identity result (99.7%), all other values placed these seven isolates as representing a novel Elizabethkingia species. We therefore propose that these seven isolates constitute a provisional novel species and propose the name Elizabethkingia umeracha; Umeracha meaning "fine waterhole" in the Peramangk language. We respectfully acknowledge the Peramangk people as the traditional owners and custodians of the waters and lands of the Adelaide Hills where these isolates originated.

Pangenome analysis
A pangenome analysis of all available Elizabethkingia genomes (n = 148) demonstrated high genetic diversity (Fig. 3). The Elizabethkingia spp. pangenome consisted of 28,240 genes, with a core genome of only 76 genes and an accessory genome of 28,164 genes (443 soft-core, 6057 shell and 21,664 cloud genes).
A pairwise genome distance MDS plot of Elizabethkingia genomes (Fig. 4) demonstrated tight clustering of all E. anophelis isolates, while E. meningoseptica isolates were the most distinct, both regarding other species and also between the E. meningoseptica isolates. The remaining Elizabethkingia species formed a more diffuse cluster with no clear distinction between human and environmental isolates and with E. umeracha sp. nov. isolates situated at the peripheries (Fig. 4, pink  triangles).
The gene presence/absence matrix generated by Roary (Supplementary Data 6) was fed into Scoary to calculate any differentiating genes present in E. umeracha sp. nov. isolates. A total of 1886 genes were only identified in these seven isolates (100% specificity, 100% sensitivity; Supplementary Data 7). More than half of these genes (n = 1110; 58.8%) encode hypothetical proteins however of the remaining genes, 537 were fed into STRING which identified several functional enrichments, the highest scoring being tryptophan biosynthesis (1.04 strength) and molybdenum cofactor biosynthesis (1.04 strength) (full analysis available in Supplementary Data 8).

Beta-lactamase resistance genes
All 94 Elizabethkingia isolates from this Australian aquatic collection carried bla B (subclass B1) and bla GOB (subclass B3) genes encoding resistance to carbapenems and a bla CME gene encoding resistance to cephalosporins.
MUSCLE alignments with all available reference sequences of bla B and bla GOB were generated to compare species and allele distributions (Fig. 6). All Elizabethkingia species in our analysis carried bla GOB , with an However, none of these deletions are expected to alter gene reading frame given they appear in multiples of three nucleotides. For the bla B gene distribution ( Fig. 6; right side tree), we saw four to five primary clades in the tree structure which were generally grouped by species.
The ten E. miricola isolates of clade 4 uniquely carried bla GOB-25-like genes and all carried bla B6-like genes. These bla B6-like genes were also identified in three Australian clinical E. miricola isolates (EkQ1, EkQ10 and EkQ13). The remaining 61 E. miricola isolates from this study uniquely carried bla GOB19-like and bla B26-like genes. Notably, E. umaracha sp. nov. isolates carried novel alleles of both metallo-β-lactamase genes.
Chromosomal extended-spectrum β-lactamase bla CME has two types known: bla CME-1 and bla CME-2. bla CME-1 appeared to be the closest allele related to the Australian aquatic environmental E. anophelis (ER-QUARD-EK_14, 22, 56) from wetland and dam samples. bla CME-2 was present at very high levels of variation from the aquatic environment isolates. Interestingly, the two distinct E. miricola clades appear to possess each a novel bla CME allele, and a third novel allele appears in the E. umeracha sp. nov. (Supplementary Data 11).

Other ARGs
No other ARGs were detected in any of the Australian aquatic environment isolates. We also searched for the known mutations in gyrA (Ser83Ile or Ser83Arg) that encode resistance to ciprofloxacin and levofloxacin (fluoroquinolone), however none were detected.

AMR phenotypic analysis: MIC testing
Ten representative isolates of E. anophelis, E. miricola and E. umeracha sp. nov., harboring unique combinations of bla B , bla GOB and bla CME genes from each Elizabethkingia clade were tested for MIC against 38 clinically relevant antimicrobials (Table 4). To date, Elizabethkingia species lack their own defined breakpoint, so they have been interpreted by using EUCAST non-species and NCSI non-Enterobacteriaceae PK-PD breakpoints (Supplementary data 12). All isolates tested showed remarkable resistance to carbapenems, cephalosporins, penicillins including carboxypenicillin and monobactam. Regarding different bla GOB, bla B and bla CME combinations, some differences in resistance profiles were noted (Table 4), including piperacillin/tazobactam resistance in only bla B-26like /bla GOB-19-like /bla CME-variant E. miricola isolates, and cefepime resistance in two E. anophelis isolates (one bla B-1-like /bla GOB-20 /bla CME-1 and one bla B-1-like /bla GOB-variant /bla CME-1 ).
Isolates were also resistant to antibiotic classes other than carbapenem and ESBLs, including aminoglycosides and glycylcycline (Table 4). One isolate, E. anophelis ER-QUAD-EK_14 was highly resistant to chloramphenicol and trimethoprim/sulfamethoxazole. Azithromycin and rifampicin have no corresponding breakpoint in EUCAST or CLSI, however, tested isolates had a very low MIC (up to the lowest range tested), suggesting a potential sensitive profile. Vancomycin and teicoplanin also lack a breakpoint but their MICs were very high, indicating non-susceptibility. Likewise for the glycopeptides and colistin, which showed higher MIC than the top of concentration tested, suggesting potential resistance of Elizabethkingia against these antibiotics.

Mobile genetic elements characterization: ICEs, plasmids and phages
Integrative conjugative elements were identified in 67 / 94 (71.3%) of the aquatic environments Elizabethkingia spp. by comparing to Elizabethkingia ICE sequences publicly available in Genbank, comprising three types of ICEs (Xu et al., 2019): ICEEaI from strain CSID3015183678, ICEEaII from strain NUHP1 and ICEEaIII from strain R26. The alignments demonstrated imperfect matches to the reference sequences, however all matches were closest to the type III ICE from E. anophelis strain R26. Two plasmids have been described from Elizabethkingia species so far (Accessions CP016375.1 and CM003640.1), however neither were detected in this study. In our analysis, both plasmid sequences were aligned to the aquatic Elizabethkingia isolates from this collection, as well clinical Elizabethkingia isolates from Australia sourced from Genbank. Thirteen (81%) of the environmental E. anophelis in clade 2 showed lowquality matches to plasmid CP016375.1 (average 8% coverage at 90% identity) (Supplementary Data 13). Alignments to the second reference plasmid from E. miricola strain EM_CHUV (CM003640.1) revealed little to no homology.   6. Elizabethkiniga bla GOB and bla B alleles. Phylogenetic trees of all Elizabethkingia bla B and bla GOB alleles. Left side is the tree of bla GOB alleles and right side is the tree of bla B alleles. Labels are colored in red for E. miricola, green for E. anophelis, blue for E. meningoceptica and orange for E. umeracha sp. nov. Connecting space between the trees links sequences from the same isolate. Available allele numbers are presented as colored strips.  We performed a phage analysis on 94 Elizabethkingia genomes via Phaster, with hits detected for each isolate (data not shown). While most hits were recorded as 'incomplete' due to scaffolding of the WGS, fourteen isolates were detected with complete phages (Supplementary Data 14). From these fourteen isolates, coming from multiple species, we selected four ranges of hit length to group the phages: 48.7 kb, 21.4-27.1 kb, 10.4 − 18.8 kb, and 4.9-9.7 kb. Alignments were generated to identify similar phages amongst the different size ranges, with only the largest (48.7 kb) appearing conserved in multiple isolates. Initially, this phage was detected in two E. miricola, one from a wetland and another from a river, with alignments demonstrating different chromosomal locations. By aligning the large phage against the Australian Elizabethkingia sequences, we identified 17 environment and clinical Elizabethkingia isolates to carry it: two E. anophelis (one from a wetland and another from a human bronchial alveolar lavage), nine E. miricola (seven from wetlands and two from a river), one E. bruuniana (human blood) and five E. umeracha sp. nov. (dam). While short-read assembly data prevents detailed comparisons of these phages, it is clear the Elizabethkingia isolates share mobile DNA.

Discussion
Elizabethkingia spp. are emerging pathogens and the only known organisms with multiple chromosomal metallo-β-lactamase genes, offering inherent resistance to carbapenems (Hu et al., 2020). Elizabethkingia species are considered environmental bacteria, with water bodies serving as environmental reservoirs. Contaminated water is implicated in Elizabethkingia spp. transmission pathways (Booth, 2014), yet with the exception of insects (Kämpfer et al., 2011), frogs (Hu et al., 2017;Lei et al., 2019), reptiles (Jiang et al., 2017) and spacecraft (Li et al., 2003), most studies on Elizabethkingia genus have focused on clinical isolates and isolates taken from hospital environments, leaving Elizabethkingia species dwelling in aquatic environments unexplored.
Here we characterized WGS of 94 Elizabethkingia derived from dam, river, and wetland samples from South Australia, thereby providing the first study of Elizabethkingia from diverse aquatic environments. Furthermore, we provide comparative genomics analyses of these environmental isolates with clinical Elizabethkingia isolates originating from Australia and worldwide.
Correctly identifying Elizabethkingia species is paramount as not only is the literature on Elizabethkingia spp. convoluted due to various nomenclature changes, but standard commercial microbial identification systems, such as biochemical tests and mass spectrometry (MS) using standard databases, cannot currently differentiate E. anophelis, E. bruuniana, E. ursingii or E. occulta and these are often misidentified as either E. meningoseptica or E. miricola Burnard et al., 2020;Lin et al., 2019b;Snesrud et al., 2019). Consistently, our initial MALDI-TOF MS results misidentified 11 isolates as E. meningoseptica. Interestingly, all speciation methods we used, including MALDI-TOF MS, Kraken2 and SpeciesFinder, gave conflicting results. The issues regarding MS and Kraken2 misidentifications likely arose from using standard databases (Lin et al., 2019a) while SpeciesFinder uses 16S rRNA gene sequences, which are known to be limited for taxonomic purposes  with studies demonstrating less than 30% accuracy for aerobic bacteria to the species level (Teng et al., 2011). These data suggest that in lieu of WGS, future Elizabethkingia spp. studies should be cautious in using 16S rRNA for speciation and ensure any utilized MS databases include all Elizabethkingia species.
The fact that E. anophelis was misidentified here and elsewhere Han et al., 2017;McTaggart et al., 2019;Snesrud et al., 2019) as E. meningosepticum by conventional clinical methods has led to speculation that E. anophelis is not only underrepresented but may actually be the primary species to cause disease in humans. This hypothesis is strengthened by recent reports of life-threatening E. anophelis infections in Asia, Australia, and the USA (Burnard et al., 2020;Lin et al., 2019a). The 16 E. anophelis isolates identified here in dam and wetland samples formed two clades which differed by approximately 807 SNPs. The single wetland isolate and two dam isolates were found to differ by only ~36 SNPs to three E. anophelis isolates originating from sepsis patients in Queensland, Australia, and ~42 SNPs to two E. anophelis isolates derived from sinks located in a Queensland hospital (Burnard et al., 2020). Screening for putative virulence factors did not identify any specific to the three environmental isolates and the three sepsis isolates, however several putative virulence factors were found to be unique to E. anophelis in general, including homologs of lipopolysaccharide biosynthesis proteins and serum resistance protein OmpA.
The majority of Australian aquatic environment isolates from this study were E. miricola. E. miricola is known to cause sepsis, oral and urinary tract infections (Green et al., 2008;Lin et al., 2019a;Zdziarski et al., 2017) though reports of infection are less frequent than that for E. meningoseptica and E. anophelis. Environmental E. miricola from this study formed two distinct clades. While a clade of ten isolates were phylogenetically positioned next to three Australian clinical E. miricola isolates taken from sputum samples (Burnard et al., 2020), the average SNP difference between these environmental and clinical isolates was 21,539. Nevertheless, environmental E. miricola shared putative virulence factors with clinical strains including homologs of capsule protein, Cps4I (Mirza et al., 2018), haemolytic toxin SmcL (González-Zorn et al., 2000), and iron acquisition protein YbtP (Fetherston et al., 1999). The remaining 61 E. miricola isolates formed a clonal clade with an average difference of 7 SNPs. Despite the highly clonal nature of these isolates, they originated from two wetland sites approximately 10 km apart. Identifying potential transmission pathways, such as through wildlife, will be critical for future epidemiology.
The remaining seven Elizabethkingia isolates originated from dam samples and were phylogenetically placed proximal to E. bruuniana strains, though SNP analyses demonstrated the two branches to differ by ~124,216 SNPs. This considerable difference prompted an investigation as to whether these seven isolates represent a novel Elizabethkingia species. WGS-based ANI and in silico DDH analyses are known to be robust speciation methods that have proven more accurate than even the gold-standard conventional DDH for bacterial species delineation (Konstantinidis and Tiedje, 2005;Lin et al., 2019a;Meier-Kolthoff et al., 2013;Richter and Rosselló-Móra, 2009;Varghese et al., 2015). Here, both analyses clearly indicated that these isolates were a distinct species. The rpoB gene is also used in speciation and possesses a higher resolution for delineation than the 16S rRNA gene (Adékambi et al., 2009). Here, all seven isolates fell under the 97.7% similarity cut-off to be classified as E. bruuniana. Furthermore, we identified 1886 genes unique to these isolates and where functional assignment was possible, these genes were mainly involved in cellular and metabolic processes. Together these data indicate a closely related but distinct species to E. bruuniana and we propose the name E. umeracha sp. nov. Though these isolates were shown here to share 62 putative virulence factors with E. anophelis and E. miricola isolates, future studies are required to determine whether this new species is pathogenic.
Metallo-β-lactamases (MBLs) are a worldwide concern as they confer resistance against carbapenems and almost all β-lactams (Chang et al., 2019), making pathogens carrying MBLs very difficult to treat. Elizabethkingia are currently the only known organisms to carry two chromosomal MBLs (bla B and bla GOB ) and additionally carry a chromosomal bla CME gene conferring resistance to cephalosporins (González and Vila, 2012). Here we found that the three species of Elizabethkingia residing in aquatic environments all carried multiple known alleles, as well as novel variants of bla B , bla GOB and bla CME genes. Regardless of allelic combinations, MIC testing demonstrated all were highly resistant to carbapenems, penicillins and monobactams. All tested isolates were also resistant to cephalosporins, though one E. miricola isolate carrying a bla CME variant was susceptible to cefepime. Intravenous vancomycin has been cited as the favorable treatment option for Elizabethkingia infections (Jean et al., 2017). While most tested isolates had MICs of vancomycin at 4 or 8 µg/mL, one E. umeracha sp. nov. isolate had > 32 µg/mL, suggesting that a conventional dose of vancomycin would not be effective (Chang et al., 2019). In addition to carbapenems and β-lactams, several isolates were found resistant to aminoglycosides and one E. anophelis isolate was highly resistant to trimethoprim/sulfamethoxazole and chloramphenicol. However, no additional AMR genes were identified, suggesting the presence of novel AMR mechanisms.
Environmental bacteria often harbor important AMR genes that are later captured and disseminated by more common human pathogens. For example, the ESBL gene, bla CTX− M , now endemic among Enterobacteriaceae, likely originated from Kluyvera ascorbate, a soil bacteria (Humeniuk et al., 2002). The acquisition and spread of resistance, as well as virulence genes, is facilitated by MGEs. ICEs are MGEs that integrate into a host chromosome and can bestow new phenotypes (Xu et al., 2019). In the 2015-2016 E. anophelis outbreak in the USA, which led to 66 confirmed cases of sepsis and 19 deaths, an ICE was identified in all of the outbreak clones. This ICE interrupted the mutY gene which led to a hypermutator phenotype (Perrin et al., 2017). ICEs were identified in 71% of 94 environmental Elizabethkingia isolates, however these bore no similarity to the ICE associated with the USA outbreak and only ~50% similarity to the type III ICE found in E. anophelis strain R26 (Kämpfer et al., 2011), the first E. anophelis strain isolated. As such we have extended knowledge on the types of ICE that are found in Elizabethkingia spp..
In conclusion, we presented the first WGS analysis of Elizabethkingia species found in aquatic environments and discovered that they carry diverse bla B bla GOB and bla CME genes and are highly resistant to carbapenems, cephalosporins, monobactams and other beta-lactams. Some isolates were also resistant to additional antibiotic classes suggesting the presence of yet undiscovered AMR mechanisms. We uncovered environmental E. anophelis isolates that were very closely related to sepsiscausing clinical strains, thus identifying water bodies as an important reservoir for pathogenic Elizabethkingia spp. and highlighting the potential for cross-habitat movement. Finally, we discovered a proposed novel species, E. umeracha sp. nov., representatives of which appear resistant to vancomycin and carry novel metallo β-lactamase and extended-spectrum serine β-lactamase gene alelles.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.crmicr.2021.100083.