The Evolutionary Trend and Genomic Features of an Emerging Lineage of Elizabethkingia anophelis Strains in Taiwan

ABSTRACT The incidence of Elizabethkingia anophelis bacteremia increased significantly in a tertiary hospital, Changhua Christian Hospital (CCH) since 2013. The infection density was 1.3 and 8.1 cases per 100,000 patient-days between 2005 and 2012 and 2013 and 2020, respectively (P < 0.05). During an outbreak investigation, a specific lineage of E. anophelis strains was identified by the pulsed-field gel electrophoresis analysis. To evaluate the evolution of the specific E. anophelis lineage, whole-genome sequencing was performed, and unique genomic features (GRs) were determined by comparative genomic analysis. The specific E. anophelis lineage was novel compared to worldwide strains ever reported by cg-MLST phylogenic and whole-genome comparative analysis. Multiplex PCR using primers designed from unique GRs were performed for prevalence screening among isolates from the CCH and nationwide isolates from the Taiwan surveillance of Antimicrobial Resistance (TSAR) Program. The proportion of the specific E. anophelis lineage increased from 7.9% (3/38) during 2005-2012 to 89.2% (223/250) during 2013-2020 (P < 0.05). Although E. anophelis usually confers resistance to multiple antibiotics with limited therapeutic options, the E. anophelis strains in the specific lineage had higher ciprofloxacin resistance (100% [226/226] versus 27.4% [17/62], P < 0.05) and was associated with a higher 14-day mortality rates (33.2% [37/226] versus 16.1% [10/62], P < 0.05) than other strains at CCH. A similarly increasing trend was also found in the national TSAR program during 2002-2018 (p for trend <0.05). We concluded that a novel lineage of E. anophelis strains has emerged dominantly in Taiwan. The genomic features are important for further investigations of epidemiology, resistance, virulence, and appropriate treatment. IMPORTANCE Elizabethkingia anophelis is an emerging multidrug resistant pathogen caused several global outbreaks recently. E. anophelis was frequently misidentified as E. meningoseptica in the past by conventional culture methods; therefore, the prevalence was often underestimated. Through revised identification, an increasing trend of E. anophelis infection was noted in a tertiary hospital and a dominant lineage of strains was recognized by genotyping. To our best knowledge, the dominant lineage of E. anophelis is novel in comparison to other worldwide strains by whole-genome comparative analysis and several unique genomic regions were found. The whole-genome sequencing data also demonstrated multiple putative virulence factors and genes associated with multidrug resistance. In our study, we identified a specially evolved E. anophelis in Taiwan with increasing nationwide dominance. This study will assist in further epidemiology surveillance and developing corresponsive infection control policies to restrain it potential of global dissemination.

IMPORTANCE Elizabethkingia anophelis is an emerging multidrug resistant pathogen caused several global outbreaks recently. E. anophelis was frequently misidentified as E. meningoseptica in the past by conventional culture methods; therefore, the prevalence was often underestimated. Through revised identification, an increasing trend of E. anophelis infection was noted in a tertiary hospital and a dominant lineage of strains was recognized by genotyping. To our best knowledge, the dominant lineage of E. anophelis is novel in comparison to other worldwide strains by whole-genome comparative analysis and several unique genomic regions were found. The wholegenome sequencing data also demonstrated multiple putative virulence factors and genes associated with multidrug resistance. In our study, we identified a specially evolved E. anophelis in Taiwan with increasing nationwide dominance. This study will assist in further epidemiology surveillance and developing corresponsive infection control policies to restrain it potential of global dissemination.

RESULTS
The trends of incidence and antimicrobial resistance of E. anophelis infections. From 2005 to 2020, a total of 334 isolates of Elizabethkingia species were retrieved at CCH. After reidentification, E. anophelis was most dominant among all Elizabethkingia spp. (86.2%, 288/334), followed by E. meningoseptica (8.7%, 29/334), E. miricola (4.8%, 16/334), and Elizabethkingia spp. (0.3%, 1/334). A significantly increasing trend of Elizabethkingia species infection was noted with increases of the case number of E. anophelis more obviously than other two species since 2013 (Fig. 1). For E. anophelis, the infection density increased from 1.3 to 8.1 cases per 100,000 patient-days between 2005 and 2012 and 2013-2020, respectively (P , 0.05). The case number of E. anophelis bacteremia peaked in 2017 with 65 cases being identified in a year. The first intrahospital outbreak was detected in 2015 in a respiratory care center with 4 patients involved. A dominant lineage of E. anophelis strains had been identified according to PFGE analysis during the early outbreak investigations (20). Intensive infection control strategies were implemented since 2016 and the case number had been reduced gradually. Only 6 cases of E. anophelis bacteremia was recorded through the whole year in 2019. However, the incidence rebounded during the COVID-19 pandemics in 2020.
The genomic features of the dominant lineage of E. anophelis strains. Twentytwo E. anophelis strains included in previous PFGE studies were selected for complete genomic sequencing and comparative genomic analysis (20). The isolates, including 18 strains belonging to the specific lineage and 4 other E. anophelis strains were selected based on their pulsotypes. The sequencing data resulted in a single circular genome for each strain with the length ranging from 3,982,150 to 4,084,052 bp and the mean GC content ranging from 35.7% to 35.9%. The cg-MLST analysis was performed on 41 completed genomes and 45 draft genomes, including our 22 isolates and other worldwide strains. The details of all genomes used in cg-MLST analysis were summarized in Table S1. From about 4,000 genes in each of the genomes, 1,779 genes were selected as core genes to generate core genome scheme and multilocus alignment. In the cg-MLST phylogenetic tree, all 18 strains in the dominant lineage and one other E. anophelis isolate (strain 277-17) were within the same clade. These 19 isolates were very closely related to each other and formed a distinct clade which was separated from other isolates with a long branch (Fig. 2). Three of our E. anophelis strains (313-22, 2-8 and 2-14) were located outside the major clade.
For further comparative genomic analysis, we used E. anophelis 296-96 as a representative strain as it was the earliest strain in the dominant lineage. The results are illustrated with the use of the BLAST Ring Image Generator (BRIG) in Fig. 3. In comparison with strain 296-96, several genomic regions (GRs) ranging from 35-kb to 90-kb were found unique among strains in the dominant lineage. Those GRs were numbered from GR1 to GR7. All 18 strains in the specific lineage had GR1 to GR6, but GR7 was only present in strains 2-62, 344-1 and 296-96. The annotated genes in the seven GRs are listed in Table 1. The seven GRs were searched among all E. anophelis strains with available genome data in NCBI database to date by the Basic Local Alignment Search Tool (BLAST). The BLAST results disclosed no other isolates had similar GRs to our novel E. anophelis strains. The only exception was E. anophelis SUE, an isolate recently recovered from another hospital in central Taiwan (21). The genome of E. anophelis SUE had all six GRs and its genome was identical to those of our strains in the dominant lineage (Fig. S1). Analysis of the E. anophelis 296-96 genome using Islandviewer4, which is empowered by integrated genomic island predictors, have shown that 5 (GR1, 3, 4, 6, and 7) out of the 7 GRs are likely genomic islands. GR1 and GR3 were predicted positive, respectively, by SIGI-HMM and IslandPath-DIMOB. GR4, GR6 and GR7 were positive by both predictors (Table S2). Annotation using PHAST shows that a 47-kb region (4,081,340-4,127,081) located within GR7 is phage-related. Integrative and conjugative elements (ICEs) are relatively common among E. anophelis genomes and the Wisconsin E. anophelis strain had it disrupted mutY gene due to insertion of an ICEEaI. Therefore, we surveyed the presence of ICEs among the genomes of our isolates. The characteristic genes associated with ICEEaIII and ICEEaI were found in GR1 and GR6, respectively. Despite an ICEEaI was found in our specific lineage, the locus of insertion was different from that of the Wisconsin strain, and the mutY gene was intact in our strains. Besides, the cargo genes carried by the ICEEaI were also very different in compared to those in the Wisconsin strain (Fig. S2). The presence and insertion locations of different ICEs of our strains and other outbreak representative strains are illustrated in Fig. S1.
Primer design and prevalence screening for outbreak strain. The primers designed for multiplex PCR are listed in Table 2 including three primer pairs were designed specifically for target sequences within GR3, GR4 and GR5. We omitted GR1 and GR6 because these are likely ICEs, which are potentially mobile and not ideal targets. GR2 was not used because it was also present in strain 277-17, a strain not belonging to the dominant lineage by PFGE analysis. To ensure the specificity of the three primer pairs, we had avoided sequences from any mobilized elements such as phages, transposons, integrons and repetitive sequences. In addition, the primers were also checked using NCBI database and their specificity were confirmed. The result of a rapid screen by multiplex PCR revealed that the E. anophelis strain with genomic features as the dominant lineage was first detected in 2011 and the case number of those specific strains increased significantly from 2011 to 2020 (Fig. 4). We also investigated 90 isolates from multiple hospitals in the TSAR program between 2002 and 2018 (18). The earliest isolate belonging to the dominant lineage was identified in 2008 from a hospital in northern Taiwan. The ratio of strains in the specific E. anophelis lineage also increased from 6.25% (3/48) during 2002-2012 to 33.3% (14/42) during 2014-2018; of note, it was 50% (n = 7/14) in 2018 (Fig. S3).
Antimicrobial susceptibility and the clinical features between patients infected by different E. anophelis strains. The antimicrobial susceptibility between strains belonging to the dominant lineage and others were compared (Table 3). All E. anophelis strains were resistant to most of the antimicrobial agents tested. The susceptible rates of trimethoprim-sulfamethoxazole were about 50-60% without statistically significant difference noted between two groups. The most prominent difference was the susceptibility to ciprofloxacin; all strains in the dominant lineage were resistant to ciprofloxacin (100%, 226/226) which was higher than other strains (27.4%, 17/62). All our 18 sequenced strains in the dominant lineage had DNA gyrase subunit A (gyrA) mutations at codons 83, Ser83Ile, that was associated with ciprofloxacin resistance. Instead, 3 other E. anophelis strains did not have the gyrA mutation (strain 313-22, 2-6 and 2-14) and another strain 277-17 carried a single mutation, Ser83Arg. All 22 sequenced E. anophelis isolates carried multiple b-lactamases genes (bla CME-1 and bla GOB ) which were associated extensive drug resistance to b-lactams. The strain 277-17 also carried genes conferring resistance to macrolides (ere(D), erm(F)), sulfonamide (sul2) and tetracycline (tet[X]) but those genes were not found among 18 sequenced strains in the dominant lineage. The MIC of minocycline determined by E-tests ranged from 0.19 to 0.5 mg/ml among 18 E. anophelis strains in our dominant lineage, which were considered susceptible (5 4 mg/ ml) by using the breakpoints for other non-Enterobacterales according to CLSI M100-S30 criteria.
For 288 patients with bacteremia due to E. anophelis, 226 (78.5%) were infected by E. anophelis strains belonging to the dominant lineage and 62 (21.5%) were by other strains ( Table 4). The average age and sex distribution were similar between two groups. Strains in the dominant lineage had higher probability to be isolated from patients admitted to respiratory care centers or respiratory wards; therefore, more patients in this group required mechanical ventilators supports (78.3% [177/226] were identified different between outbreak and non-outbreak strains. GR1 and GR6 were type III and type I integrative and conjugative elements (ICEs), respectively. versus 46.8%, [29/62], P , 0.001). Most patients were treated with trimethoprim-sulfamethoxazole. Fluoroquinolones were less used for patients infected by strains in the dominant lineage because of high resistance rates; instead, more patients received minocycline as alternative treatment (Fig. 4). Besides, glycopeptide (vancomycin) and penicillin (piperacillin) were more commonly used for patients infected by other strains in early period, which was probably attributable to former susceptibility reports determined by disk methods showing inhibitory zones for vancomycin and piperacillin before 2010. Finally, the mortality rates were higher among patients infected by strains in the dominant lineage than those by other strains both in 14-day mortality rates (33.2%, 75/226 versus 16.1%, 10/62, P = 0.009) and in-hospital mortality rates (50.9%, 115/226 versus 32.2%, 20/62, P = 0.014).

DISCUSSION
In this study, we identified a specific lineage of E. anophelis strains that evolved and emerged in recent decades in Taiwan. Through whole-genome sequencing and comparative genomic analysis, the specific E. anophelis lineage were proved different from any currently known isolates around the world and associated with extensive antimicrobial resistance and high mortality rates.  A recent study explored the relatedness of 22 indigenous E. anophelis in Australia over a period of 16 years. Some E. anophelis isolates in Australia are genetically related to the strains from the Unites States, England, Singapore and Taiwan, which implies international spread of E. anophelis isolates (15). From 2015 to 2016, a unique E. anophelis strain was reported in a large-scale outbreak in WI, USA (7,10). Similarly, a specific lineage of E. anophelis strains was identified in our study with unique genomic features (20). Although cg-MLST genotyping discriminated our E. anophelis strains from the others well, it has not escaped our notice that, strain 277-17, which does not belong to the dominant lineage by PFGE analysis was clustered together with the 18 dominant lineage stains in the cg-MLST tree (22). It is the comparative genomics analysis on the presence of the GRs   strains in the dominant lineage were lack of tet(X) gene that associated with tetracycline resistance. The MICs of 18 strains were all less than 4 mg/ml. The result was similar to those of many recent studies (24)(25)(26). A study in Shanghai, China evaluated antimicrobial susceptibility of 52 Elizabethkingia isolates (including 35 E. anophelis, 14 E. meningoseptica, and 2 E. miricola) and minocycline was more active than doxycycline and tigecycline (# 4 mg/ml, 100% versus 96.2% and 78.8%, respectively) (24). Another study in Taiwan also demonstrated that 100% (90/90), 96.7% (87/90), and 52.2% (47/ 90) of E. anophelis isolates were susceptible to minocycline, doxycycline and tigecycline, respectively (18). Though no current breakpoint value is available, the clinical use of tetracycline, especially for minocycline and doxycycline, to treat E. anophelis infections warrants more investigations. Several gene clusters related to membrane-associated transporter functions were identified in the GRs specific to the specific lineage of E. anophelis strains; these included genes for an MFS-type transporter in GR2, genes for an RND-type efflux system in GR4, gene clusters for two ABC-type transporter systems in GR6, genes for outer membrane receptors and several genes for porins. The MFS, RND, and ABC-type transporters have been long regarded as members of the bacterial efflux pump systems that may play a role in facilitating the efflux transport of antimicrobials, ions, or bioactive molecules across the membrane (27,28). In a recent review, the bacterial efflux pump is not only associated with antibiotic resistance but also play roles in quorum sensing, biofilm formation, pathogenicity and virulence (29,30). It is plausible that these extra transporter genes in the GRs may be a benefit to the dominant of the specific lineage of E. anophelis strains.
There are several limitations in this study. First, molecular epidemiologic studies could be extended to investigate the local and global trends of E. anophelis infections related to the specific lineage with time. Second, while some antimicrobial agents such as vancomycin or rifampin were reported to be active in vitro against some Elizabethkingia species (18, 31), we did not assess the impact of antimicrobial therapy chosen based on the available antimicrobial susceptibility on the outcome of E. anophelis infections.
In conclusion, our study revealed an increasing trend of infection due to a specific lineage of E. anophelis strains. Those specific strains were associated with extensive antimicrobial resistance and high mortality rates. The genomic features identified would be helpful for further evaluation of the evolving epidemiology, antimicrobial resistance, virulence and modes of transmission.

MATERIALS AND METHODS
Setting and bacterial isolates. All nonduplicated Elizabethkingia species isolated from patients between 2005 and 2020 were retrieved retrospectively from the biobank at Changhua Christian Hospital (CCH). All isolates were reidentified by the Matrix Assisted Laser Desorption Ionization Time-of-Flight (MALDI-TOF) technology, VitekMS with Knowledge Base Version 3.0 (bioMérieux SA, Marcy L' Etoile, France) as previously reported (17,20). The antimicrobial susceptibility was determined using Vitek 2 system and susceptibility testing cards (bioMérieux, Lyon, France). The MIC of minocycline was determined by Etest (bioMérieux, Marcy l'Etoile, France) according to the manufacturer's recommendations (32). The breakpoints for antimicrobial susceptibility was interpreted using criteria for other non-Enterobacterales according to CLSI M100-S30 (33). The breakpoint used for colistin was adapted from the criteria for Pseudomonas aeruginosa that a MIC $ 4 mg/ml is categorized as resistance in the CLSI M100-S30 (33). A standardized case record form was used to collect information on patients with E. anophelis bacteremia.
The E. anophelis isolates collected in the Taiwan Surveillance of Antimicrobial Resistance (TSAR) program from 2002 (TSARIII) to 2018 (TSARIX) were included in trend analysis. The TSAR program is a biannual nationwide bacteria repository program in Taiwan since 1998 (18).
Whole-genome sequencing, core-genome multilocus sequence typing, and comparative genomic analysis of E. anophelis strains. Whole-genome sequencing were performed using Illumina iSeq100 for short-read sequencing and Nanopore MinION platforms for long-read sequencing as reported previously (34). Briefly, the libraries of short-read sequencing were constructed using the Nextera DNA CD Indexes (Illumina) and analyzed by Agilent 2100 Bioanalyzer that generated 150 bp paired end reads, followed by adaptor trimming. For long-read sequencing, libraries were constructed using Rapid Barcoding kit SQK-RBK004 (Oxford Nanopore Technologies) and sequenced on a Nanopore MinION following standard protocols (34). Hybrid assembly of the Illumina short paired-end reads and Nanopore long reads were performed using Unicycler version 0.4.8 to result in a single circular chromosome of all isolates (4). The correctness of the assembly results was rechecked by using another assembler, Canu 1.8, with the aid of CLC genomic workbench 11 (Qiagen). The genomes were annotated by the NCBI Prokaryotic Genome Annotation Pipeline (PGAP). The core-genome multilocus sequence typing (cg-MLST) phylogenetic tree analysis was conducted using Harvest suite, Parsnp v1.2 module with default settings (35). The phylogenetic result was visualized using CLC genomic workbench v11.0.11 (Qiagen).
Comparative genomic analyses were performed using blast, MAUVE v1.0 and BRIG 0.95 to determine strain-specific genomic regions (36,37). Antimicrobial resistance gene were predicted using Resfinder v3.2 (38). Three types of integrative and conjugative elements (ICEs) were determined according to the architecture of the conjugation module and the phylogeny of the relaxase, coupling protein, TraG, and TraJ protein sequences as reported previously (35). Islandviewer4 with two genomic island prediction tools (SIGI-HMM and IslandPath-DIMOB) was used to predict genomic islands (39). Phage genomes were analyzed by PHAST (40). With reference to the result of comparative genomic analysis, we designed primer pairs in CLC genomic workbench 11 for rapid epidemiological survey. Oligomers were synthesized by Tri-I Biotech Inc. (Taichung city, Taiwan).
Ethics and experimental biosafety statements. This study was approved by the Institutional Review Board of Changhua Christian Hospital (registration no. 180502). Informed consent was waived. The experiments were also approved by the Institutional Biosafety Committee of Changhua Christian Hospital.
Statistical analysis. MedCalc software (MedCalc Software Ltd., Los Angeles, CA, USA) was used for statistical analyses. Chi-square or Fisher's exact tests were used to compare categorical variables. Noncategorical variables were compared by Student's t test for data with normal distribution or by the Mann-Whitney U-test for data without normal distribution. Trend analysis was performed by Cochran-Armitage test. All tests were two-tailed and a P value of ,0.05 was considered significant.
Data availability.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.5 MB.