microorganisms Article Comparative Genomic Analyses of Virulence and Antimicrobial Resistance in Citrobacter werkmanii, an Emerging Opportunistic Pathogen José R. Aguirre-Sánchez 1 , Beatriz Quiñones 2 , José A. Ortiz-Muñoz 3 , Rogelio Prieto-Alvarado 3 , Inés F. Vega-López 3 , Jaime Martínez-Urtaza 4 , Bertram G. Lee 2 and Cristóbal Chaidez 1, * 1 Laboratorio Nacional para la Investigación en Inocuidad Alimentaria, Centro de Investigación en Alimentación y Desarrollo A.C. (CIAD), Coordinación Regional Culiacán, Culiacan 80110, Mexico; jose.aguirre.dc18@estudiantes.ciad.mx 2 Produce Safety and Microbiology Research Unit, Western Regional Research Center, Agricultural Research Service, U.S. Department of Agriculture, Albany, CA 94710, USA; beatriz.quinones@usda.gov (B.Q.); bertram.lee@usda.gov (B.G.L.) 3 Parque de Innovación Tecnológica de la Universidad Autónoma de Sinaloa, Culiacan 80040, Mexico; armando_3133@outlook.com (J.A.O.-M.); rogelioprieto@uas.edu.mx (R.P.-A.); ifvega@uas.edu.mx (I.F.V.-L.) 4 Departament de Genètica i de Microbiologia, Universitat Autờnoma de Barcelona, 08193 Bellaterra, Spain; jaime.martinez.urtaza@uab.cat * Correspondence: chaqui@ciad.mx; Tel.: +52-(667)-480-6950 Abstract: Citrobacter werkmanii is an emerging and opportunistic human pathogen found in devel- oping countries and is a causative agent of wound, urinary tract, and blood infections. The present study conducted comparative genomic analyses of a C. werkmanii strain collection from diverse geographical locations and sources to identify the relevant virulence and antimicrobial resistance Citation: Aguirre-Sánchez, J.R.; genes. Pangenome analyses divided the examined C. werkmanii strains into five distinct clades; Quiñones, B.; Ortiz-Muñoz, J.A.; the subsequent classification identified genes with functional roles in carbohydrate and general Prieto-Alvarado, R.; Vega-López, I.F.; metabolism for the core genome and genes with a role in secretion, adherence, and the mobilome for Martínez-Urtaza, J.; Lee, B.G.; the shell and cloud genomes. A maximum-likelihood phylogenetic tree with a heatmap, showing the Chaidez, C. Comparative Genomic virulence and antimicrobial genes’ presence or absence, demonstrated the presence of genes with Analyses of Virulence and functional roles in secretion systems, adherence, enterobactin, and siderophore among the strains Antimicrobial Resistance in belonging to the different clades. C. werkmanii strains in clade V, predominantly from clinical sources, Citrobacter werkmanii, an Emerging Opportunistic Pathogen. harbored genes implicated in type II and type Vb secretion systems as well as multidrug resistance Microorganisms 2023, 11, 2114. to aminoglycoside, beta-lactamase, fluoroquinolone, phenicol, trimethoprim, macrolides, sulfon- https://doi.org/10.3390/ amide, and tetracycline. In summary, these comparative genomic analyses have demonstrated highly microorganisms11082114 pathogenic and multidrug-resistant genetic profiles in C. werkmanii strains, indicating a virulence potential for this commensal and opportunistic human pathogen. Academic Editors: Shengxi Chen and Fabio Zicker Keywords: Citrobacter; nosocomial infections; virulence; antimicrobial resistance; comparative Received: 14 July 2023 genomics; emerging pathogens; pangenomics; whole-genome sequencing Revised: 11 August 2023 Accepted: 13 August 2023 Published: 19 August 2023 1. Introduction The genus Citrobacter, a member of the family Enterobacteriaceae, is composed of Gram- Copyright: © 2023 by the authors. negative bacteria that are non-spore-forming bacilli and that have the ability to use citrate Licensee MDPI, Basel, Switzerland. as a carbon source [1–3]. To date, the genus has been classified into 18 species, based This article is an open access article on traditional and molecular techniques [4–6]. Citrobacter species have many reservoirs, distributed under the terms and including the human and animal gastrointestinal tracts, and these species can also be found conditions of the Creative Commons in water, soil, and food [1–3]. Several transmission routes have been proposed for this Attribution (CC BY) license (https:// bacterial pathogen, such as fecal-oral transmission, contaminated food, hospital equipment, creativecommons.org/licenses/by/ and person-to-person transmission. 4.0/). Microorganisms 2023, 11, 2114. https://doi.org/10.3390/microorganisms11082114 https://www.mdpi.com/journal/microorganisms Microorganisms 2023, 11, 2114 2 of 17 Citrobacter is considered an opportunistic pathogen and has been implicated as a causative agent of hospital settings (nosocomial) and community-acquired infections among immunocompromised patients and neonates [2,7–10]. Disease symptoms in hu- mans that are caused by Citrobacter include urinary tract infections, bloodstream infections, brain abscesses, respiratory tract infections, and neonatal infections such as meningitis and bacteremia. When compared to other Enterobacteriaceae pathogens, Citrobacter species are considered to have low virulence since human infections are uncommon in the general population [11,12]. However, evidence revealed that Citrobacter species were responsi- ble for 3–6% of all nosocomial infections attributed to the Enterobacteriaceae pathogens in surveys conducted in North America [3]. Among the Citrobacter species, Citrobacter freundii and Citrobacter koseri are the species most commonly implicated in causing the majority of opportunistic human infections [11,12]. Given that Citrobacter can persist in a host for long periods [13], Citrobacter infections can result in multidrug resistance outbreaks among neonates and immunocompromised patients who have prolonged hos- pital stays [8,10,14–17]. Recent evidence has demonstrated that other Citrobacter species, such as Citrobacter werkmanii, are also considered emerging opportunistic pathogens in developing countries [9,18–20]. Although considered to be commensal in humans and animals, C. werkmanii has previously been linked as a causative agent of wound infection, urinary tract infections, and bacteremia in humans [18,19]. More recently, the use of whole-genome sequencing has revealed that the C. werkmanii strain AK-8, isolated from a patient suffering chronic kidney disease, harbors genes with a functional role in virulence and antimicrobial resistance [20], highlighting the need for subsequent studies to further characterize the virulence potential of this commensal as an opportunistic human pathogen. To expand on the characterization of C. werkmanii as an emerging opportunistic pathogen, the present study conducted comparative genomic analyses of C. werkmanii strains recovered from a major agricultural region producing horticultural products in Mexico, as well as from various geographical locations and sources, to identify relevant virulence and antimicrobial resistance genes. This research identified key determinants implicated in highly pathogenic and multidrug- resistant profiles in C. werkmanii strains and has provided fundamental information to enable the characterization of this emerging and opportunistic pathogen in humans. 2. Materials and Methods 2.1. Bacterial Strain Isolation and Growth Conditions The examined Citrobacter strains, LANIIA-031 and LANIIA-032, were isolated from a field survey study conducted in major agricultural rivers in the state of Sinaloa in North- western Mexico [21]. Specifically, a size-exclusion ultrafiltration method was employed to efficiently concentrate large river water volumes, resulting in a final suspension of the targeted bacterial species in the retentive volume, as recently documented in [21]. The strains named LANIIA-031 and LANIIA-032 were recovered after subjecting the concen- trated river water samples to an enrichment step in tryptic soy broth (Becton Dickinson Bioxon, Mexico City, Mexico) at 37 ◦ C for 24 h, followed by subsequent growth on xylose lysine desoxycholate selective agar, as described in a previous study [21]. The strains were subsequently preserved at −80 ◦ C in a 50/50 mixture of glycerol and tryptic soy broth for further genomic characterization. 2.2. Genome Sequencing and Annotation To characterize the recovered Citrobacter strains from river water, the LANIIA-031 and LANIIA-032 strains were subjected to whole-genome sequencing. DNA extractions were performed using the DNeasy Blood and Tissue kit (QIAGEN, Mexico City, Mexico), following the manufacturer’s specifications, and the recovered genomic DNA was ini- tially assessed and quantified with a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). For performing the sequencing reactions, genomic DNA from the Citrobacter strains was quantified with a Qubit™ 2.0 Fluorometer (Invitrogen, Microorganisms 2023, 11, 2114 3 of 17 Carlsbad, CA, USA) and adjusted to a 0.2 ng/µL concentration. The genomic DNA libraries per strain were prepared using a final amount of 1 ng with the Nextera XT DNA Library Preparation Kit (Illumina Inc., San Diego, CA, USA). They were then sequenced using a MiSeq™ Reagent Kit v2 (300-cycle format) to obtain a 2 × 150 bp paired-end read output with a MiSeq™ System (Illumina, Inc.) at the Earlham Institute (Norwich Research Park, Norwich, United Kingdom), as in previous studies [22]. To initially assess the quality of the sequencing output prior to assembly, the raw sequence data were visualized and evalu- ated using the FastQC program [23]. The script wrapper Trim-Galore, version 0.6.4 [24], was employed for the removal of low-quality bases (value < 30) from the 30 -end of the reads, of adapter sequences, and of reads shorter than 80 bp in length, and the Clumpify tool, version 38.75 [25], was used to eliminate the duplicated sequencing reads. For the de novo read assembly, the A5-miseq pipeline, version 20160825 [26], was used with the re-scaffolding process to reduce the contig numbers, and the results were chosen according to the final number of contigs, N50, and genome size. To perform genome annotation, Prokka software, version 1.14.5 [27] and the RAST server [28] for the Citrobacter genus (NCBI:txid544) were used with a total of 30 genomic sequences, downloaded from the National Center for Biotechnology Information, detailing C. werkmanii strains from various geographical locations and sources. To screen the contig assemblies for the examined C. werkmanii strains, the virulence finder database and the comprehensive antibiotic resistance database [29,30] for genes related to virulence and antimicrobial resistance traits, respectively, were employed in conjunction with the ABRicate tool, version 1.0.1, by selecting the parameter cutoffs of 90% coverage and 95% nucleotide identity [31]. The ABRicate results were then repre- sented on a clustered heatmap, depicting the presence or absence of the virulence and antimicrobial resistance gene profiles, and the heatmap was constructed and edited with the Interactive Tree Of Life (iTOL) Annotation Editor, version 5 [32]. Genes encoding the proteins linked to the various bacterial secretion systems were detected by employing the program MacSyFinder, version 1, with the TXSScan model, version 1.1.1 [33,34], by using a maximal E-value of 0.001 as the statistical threshold of significance for the analysis of secretion systems. 2.3. Phylogenetic, Comparative Genomics, and Pangenome Analyses The phylogenetic analysis of the C. werkmanii strains, LANIIA-031 and LANIIA- 032, recovered from river water was performed using the concatenated sequences of the housekeeping genes, fusA (protein synthesis elongation factor-G), leuS (leucine tRNA syn- thetase), pyrG (CTP synthetase), and rpoB (β-subunit of RNA polymerase), as an established sequence-based method for the identification of species within the genus Citrobacter [35]. The sequence of recN (DNA repair) was also employed as another reliable marker for molecular species identification [36]. For the phylogenetic analyses, the genome data from the Citrobacter type strains (Table S1 in the Supplementary Materials) were also included in the phylogenetic analyses. Sequence concatenation and alignments were performed with Geneious software, version 9.1.8 (Biomatters, Auckland, New Zealand) and were im- ported into MEGA-X [37] for constructing the phylogeny trees, using the neighbor-joining method [38]. The phylogenetic tree was rooted at the midpoint, then the topology was validated by performing a bootstrap test for a total of 1000 replicates. The evolutionary dis- tances were computed using the maximum composite likelihood method, with a uniform rate for substitution [39]. Based on the results of the phylogenetic trees, the genomic sequences from closely re- lated Citrobacter species were subsequently selected for further comparison with the strains LANIIA-031 and LANIIA-032, based on documented methods for phylogenetic analysis us- ing the complete genome data [40]. Digital DNA–DNA hybridization was performed with the server located at: https://www.dsmz.de/services/online-tools/genome-to-genome- distance-calculator-ggdc (accessed on 9 December 2022), using “suggest method 2” [41]. The average nucleotide identity values between genomes and coverage were determined Microorganisms 2023, 11, 2114 4 of 17 using JSpecies, version 1.2.1 [42], and the reported average nucleotide identity values were adjusted by the genome coverage on the JSpecies web server. To initially examine genomic variability in the recovered strains of LANIIA-031 and LANIIA-032, the genome assemblies were submitted to the IslandViewer 4 server [43] to search for genomic islands. To further characterize the genomic islands, the genomes were examined for identifying phages using PHASTER [44]. Comparison of the genomes of the strains LANIIA-031 and LANIIA-032 was performed using progressiveMauve [45], and the diagram with the genomic islands was generated with Geneious software, version 9.1.8 (Biomatters). To examine the complete genetic composition, using the available sequence data of the C. werkmanii strain collection for the phylogenetic clade grouping, a pangenome was constructed using a combination of Roary, version 3.11.2 [46], and Anvi’o, version 7 [47] platforms. Phylogenetic relationships based on the core genomes of C. werkmanii were generated with the HarvestTools suite, version 1.2.2, by employing the multi-aligner, Parsnp, version 1.7.4, and visual platform, Gingr, version 1.3, to obtain a multi-FASTA alignment [48]. The RAxML program, version 8 [49], with a general time-reversible model, was used to construct a core genome phylogenetic tree using a gamma distribution and 1000 bootstrap replicates. The iTOL program, version 5 [32] was used to visualize and edit the core phylogenetic tree. Additionally, GFF files were generated with the Prokka tool, version 1.14.5; these files were used as input for generating a gene alignment and for iden- tifying the presence of core and accessory genes, based on a 90% identity with the Roary program, version 3.11.2. The pangenome was visualized with Phandango, version 1.3.0 [50], and some additional plots were generated using the script entitled “roary_plots.py”. More- over, the Anvi’o platform, version 7, was used by selecting the pangenomics workflow system (available at https://merenlab.org/2016/11/08/pangenomics-v2/ and accessed on 18 November 2021) for constructing the pangenome of the C. werkmanii strains. The results obtained using the annotated genomic databases were subsequently visualized using the command entitled “Anvi-display-pan”. The gene calls were clustered into bins, based on the following criteria: (i) core genes present in 99–100% of the genomes; (ii) shell genes present in 15–99% of the genomes; and (iii) cloud genes present in less than 15% of the genomes [46,51]. Furthermore, the data in the gene cluster summary file from the Anvi’o workflow system was used to analyze the distribution of functional clusters of orthologous gene (COG) categories in the shell and cloud pangenome in the examined C. werkmanii strains. 2.4. Statistical Analyses Statistical analyses were performed by conducting Fisher’s exact test, using the R Statis- tical Software (version 4.2.0; R Foundation for Statistical Computing, Vienna, Austria) [52]. Probability values (p-values) of lower than 0.05 were considered significant. 3. Results 3.1. Characterization of Citrobacter Strains from River Water Phylogenetic analyses of Citrobacter reference strains (Table S1 in the Supplementary Materials), based on concatenation of the MLST genes fusA, leuS, pyrG, and rpoB (Figure 1a) and also of the housekeeping recN gene (Figure 1b), revealed that the river water strains LANIIA-031 and LANIIA-032 belong to the C. werkmanii species and were found to cluster together in a separate branch with the reference strain C. werkmanii FDAARGOS_364T as the most closely related species. Based on the results of the phylogenetic trees (Figure 1), the genomic sequences of the Citrobacter river water strains, LANIIA-031 and LANIIA-032, were further examined by calculating the digital DNA–DNA hybridization, in compari- son to the genome data from closely related Citrobacter reference strains (Table S1 in the Supplementary Materials). As shown in Table 1, the recorded values of the digital DNA– DNA hybridization analysis revealed that strains LANIIA-031 and LANIIA-032 had values of >90% when compared to the reference strain, C. werkmanii FDAARGOS_364T . In contrast, the digital DNA–DNA hybridization values of strains LANIIA-031 and LANIIA-032 when reference strains (Table S1 in the Supplementary Materials). As shown in Table 1, the recorded values of the digital DNA–DNA hybridization analysis revealed that strains Microorganisms 2023, 11, 2114 LANIIA-031 and LANIIA-032 had values of > 90% when compared to the reference strain, 5 of 17 C. werkmanii FDAARGOS_364T. In contrast, the digital DNA–DNA hybridization values of strains LANIIA-031 and LANIIA-032 when compared to other species were below 70%, which is the established threshold for bacterial species delineation [5,53]. As an additional compared to other species were below 70%, which is the established threshold for bacterial in silico method for corroborating the speciation of the river water strains, LANIIA-031 species delineation [5,53]. As an additional in silico method for corroborating the speciation and LANIIA-032, the genomes were compared via the average nucleotide identity method of the river water strains, LANIIA-031 and LANIIA-032, the genomes were compared via (Table S2 in the Supplementary Materials), and the results showed values above 90% for the average nucleotide identity method (Table S2 in the Supplementary Materials), and the the LANIIA-031 results showedand LANIIA-032 values above 90% strains when for the compared LANIIA-031 andtoLANIIA-032 the reference strain, strains whenC. werkmanii FDAARGOS_364 T. When compared to other Citrobacter species, values below T compared to the reference strain, C. werkmanii FDAARGOS_364 . When compared to other 79% were observed; Citrobacter these findings species, values werewere below 79% in agreement observed;with thefindings these phylogenetic analysis were in and agreement digital DNA–DNA with the hybridization, phylogenetic analysis andindicating that the river digital DNA–DNA water strains hybridization, of LANIIA-031 indicating that the and LANIIA-032 are C. werkmanii. river water strains of LANIIA-031 and LANIIA-032 are C. werkmanii. (a) (b) Figure 1. Phylogenetic relationships of the Citrobacter strains, LANIIA-031 and LANIIA-032, Figure 1. Phylogenetic recovered relationships from agricultural of theto river water, Citrobacter strains, LANIIA-031 other Citrobacter-type strains.and LANIIA-032, recov- Maximum-likelihood ered from agricultural phylogenetic trees wereriver constructed, basedCitrobacter-type water, to other strains.fusA, on the concatenated Maximum-likelihood leuS, pyrG, and rpoBphylogenetic genes (a) and the trees recNconstructed, were gene (b), bybased usingonthethe neighbor-joining method, concatenated fusA, leuS,rooted at the pyrG, and midpoint, rpoB and genes (a) validated and the recN by performing a bootstrap test for a total of 1000 replicates. Bootstrap values (%) are indicated gene (b), by using the neighbor-joining method, rooted at the midpoint, and validated by performing at the nodes, and the scale bar represents the expected number of nucleotide substitutions. The superscript a bootstrap test for a total of 1000 replicates. Bootstrap values (%) are indicated at the nodes, and the “T” indicates a Citrobacter type strain. scale bar represents the expected number of nucleotide substitutions. The superscript “T” indicates a Citrobacter type strain. Microorganisms 2023, 11, 2114 6 of 17 Table 1. Digital DNA–DNA hybridization of the Citrobacter strains examined in the present study. C. werkmanii FDAARGOS_364 T C. werkmanii LANIIA-032 C. youngae CCUG30791 T C. braakii ATCC 51,113 T C. freundii ATCC 8090 T C. tructae SNU WT2 T C. portucalensis A60 T C. pasteurii CIP55.33T C. europaeus 97/79T C.cronae Tue2-1T Citrobacter Species 1 C. tructae SNU WT2 T 32.6 C. portucalensis A60 T 33.0 42.4 C. europaeus 97/79 T 50.3 32.9 42.9 C. braakii ATCC 51113 T 52.9 48.5 33.1 42.7 C. pasteurii CIP 55.13 T 38.8 38.5 38.4 32.9 36.2 C. youngae CCUG 30791 T 59.6 39.0 39.0 39.4 32.9 36.5 C. freundii ATCC 8090T 39.8 35.2 45.0 45.2 52.4 29.8 37.6 C. werkmanii 37.8 36.4 36.1 42.7 43.4 42.4 32.4 70.0 FDAARGOS_364 T C. werkmanii LANIIA-032 92.2 37.7 36.5 36.1 42.6 43.3 42.5 32.4 70.0 C. werkmanii LANIIA-031 92.3 92.2 37.9 36.5 36.2 42.7 43.3 42.5 32.4 70.0 1 The superscript “T” indicates a Citrobacter-type strain. 3.2. Pangenome Analyses of C. werkmanii Strains from Diverse Sources and Locations To obtain a better understanding of the virulence potential in C. werkmanii, a pangenome analysis was conducted by employing the Roary and Anvi’o workflows with an additional thirty publicly available genomic data points from C. werkmanii strains (Table 2), recov- ered from distinct geographical locations and sources. For this pangenome analysis, the resulting gene calls were clustered into bins based on: (i) core genes present in 99–100% of the genomes; (ii) shell genes present in 15–99% of the genomes; and (iii) cloud genes present in less than 15% of the genomes. A total of 11,680 gene clusters with 153,862 genes were detected to be specific for the C. werkmanii pangenome, which was classified as a core, shell, and cloud genome. As indicated by the outer colored ring in the pangenome diagram (Figure 2), the core genome in the examined C. werkmanii strains was the most abundant and consistently detected fraction (99–100% of the strains), with 3871 gene clusters and 126,303 gene calls. Moreover, 367 gene clusters and 6832 gene calls were detected for the shell genome (15–95% of the strains). Finally, the cloud genome (0–15% of the strains) accounted for 7442 gene clusters and 20,727 gene calls. The construction of a Heap’s Law chart showed that as new genomes were added to the analysis, the number of conserved genes (the core genome) decreased slightly (Figure 2, solid line), but the total number of genes increased considerably (Figure 2, dashed line). These observations were indicative of an open pangenome for C. werkmanii where the genomic content increased, with the number of additionally sequenced strains contributing to the species characterization. Microorganisms 2023, 11, 2114 7 of 17 Table 2. List of Citrobacter werkmanii strains examined in the present study and their characteristics. Clade 1 Strain 2 GenBank Accession Number Isolation Source Description Sample Source Type Date 3 Country Clade I CB00011 GCA_016505505.1 Sputum Human 2017 United States CB00078 GCA_016507625.1 Wound Human 2018 United States Clade II FDAARGOS_364 T GCA_002386385.1 Stool Human 2014 United States MGYG-HGUT-02535 GCA_902388105.1 Gut Human 2019 United States Clade III CB00044 GCA_016505055.1 Not collected Human 2017 United States CVM 45625 GCA_015942525.1 Unknown Environmental 2019 United States Clade IV BF-6 GCA_002025225.1 Industrial water Environmental 2012 China Colony242 GCA_016893825.1 Food Food 2019 Thailand Colony247 GCA_016893645.1 Food Food 2019 Thailand CVM 45620 GCA_015943365.1 Unknown Environmental 2019 United States CVM 45667 GCA_015943485.1 Unknown Environmental 2019 United States CVM 45672 GCA_015943405.1 Unknown Environmental 2020 United States LANIIA-031 JAUJUK000000000 River water Environmental 2018 Mexico LANIIA-032 JAUJUL000000000 River water Environmental 2018 Mexico RS189 GCA_015958985.1 Not collected Human 2017 United States YDC730 GCA_015958865.1 Pelvic abscess Human 2015 United States Clade V 2580 GCA_009907085.1 Urine Human 2015 Nigeria AK-8 GCA_002114305.1 Human urine Human 2014 India BOX-5 GCA_009856875.1 Hospital sink Environmental 2016 France C5.1 GCA_008364715.1 Sprouts Food 2015 Germany CRE1173 GCA_018106225.1 Pus Human 2015 Malaysia CRE1470 GCA_018106165.1 Peritoneal fluid Human 2016 Malaysia CRE346 GCA_018106145.1 Foot ulcer Human 2015 Malaysia CRE806 GCA_018106185.1 Foot ulcer Human 2014 Malaysia CRK0001 GCA_002185305.2 Blood Human 2014 United States CW_LB-887 GCA_013303045.1 Coastal water Environmental 2014 Brazil FDAARGOS_616 GCA_008693645.1 Clinical isolate Human NA United States ICR003007 GCA_004146135.1 Hospital patient Human 2017 France JS97 GCA_009821535.1 Chicken Animal NA Unknown L38 GCA_013618825.1 Chicken liver Animal 2019 Nigeria UMH18 GCA_003665555.1 Bacteremia Human 2013 United States YDC667-1 GCA_013336965. Lung tissue Human 2014 United States 1 Clade grouping was based on whole genome analysis of C. werkmanii strains, as described in Section 2, Materials and Methods. 2 The superscript “T” indicates a Citrobacter-type strain. 3 Date refers to the sample collection date. NA refers to data that are not available. Microorganisms 2023, 11, x FOR PEER REVIEW 8 of 17 Microorganisms 2023, 11, 2114 8 of 17 Figure 2. Pangenome analysis of the examined Citrobacter werkmanii strains. Pangenome analysis of Figure 2. Pangenome over thirty analysis C. werkmanii of the examined sequenced Citrobacter strains (Table 2) waswerkmanii strains. constructed usingPangenome Roary and analysis Anvi’o ofplatforms, over thirty C. werkmanii and the resulting gene calls were clustered into bins, based on: (i) coreRoary sequenced strains (Table 2) was constructed using genesand Anvi’o present in 99–100% of platforms, and thethe genomes; resulting(ii)gene shellcalls genes present were in 15–99% clustered of the into bins, genomes; based on: (i) and core (iii) genes cloud genes present in present in less than 15% of the genomes. The construction of a phylogram, based 99–100% of the genomes; (ii) shell genes present in 15–99% of the genomes; and (iii) cloud geneson gene presence or absence present in lessinthan the examined 15% of theC. werkmanii genomes. strains, The resulted construction of in the grouping a phylogram, of the based on strains into five gene presence or distinct clades. Construction of a Heap´s Law chart (top right) showed that the number of conserved absence in the examined C. werkmanii strains, resulted in the grouping of the strains into five distinct genes in the core genome decreased slightly (solid line) but the total number of genes increased clades. Construction of a Heap´s Law chart (top right) showed that the number of conserved genes in considerably (dashed line), indicative of an open pangenome for C. werkmanii. the core genome decreased slightly (solid line) but the total number of genes increased considerably (dashed line), indicative Phylogenetic of an open analysis pangenome of the for C. werkmanii. C. werkmanii strains divided the pangenome of this species into five distinct clades (Figure 2 and Table 2). The phylogram, based on gene Phylogenetic analysis of the C. werkmanii strains divided the pangenome of this species presence or absence, revealed that most examined C. werkmanii strains were clustered in into five distinct clades (Figure 2 and Table 2). The phylogram, based on gene presence two major clades, namely, clades IV and V. In particular, clade IV was composed of ten or absence, revealed that most examined C. werkmanii strains were clustered in two major strains, including the two river-derived C. werkmanii strains, LANIIA-031 and LANIIA- clades, namely, clades IV and V. In particular, clade IV was composed of ten strains, includ- 032, from this study, as well as strain BF-6 from industrial water, strains Colony242 and ing the two river-derived C. werkmanii strains, LANIIA-031 and LANIIA-032, from this Colony247 from food, and the environmental strains CVM45667, CVM45620, and study, as well as strain BF-6 from industrial water, strains Colony242 and Colony247 from CVM45672. A correlation analysis revealed that this clade was significantly associated food, and the environmental strains CVM45667, CVM45620, and CVM45672. A correlation with environmental sources (Fisher’s exact test p-value of < 0.01). Moreover, clade V, with analysis revealed that this clade was significantly associated with environmental sources the largest number of strains (Figure 2 and Table 2), a total of 16 strains, comprised C. (Fisher’s exact test p-value of < 0.01). Moreover, clade V, with the largest number of werkmanii strains that had predominantly been recovered from clinical samples, strains (Figure 2 and Table 2), a total of 16 strains, comprised C. werkmanii strains that associated with minor and severe human illnesses. The clinical strains included AK-8, had predominantly been recovered from clinical samples, associated with minor and se- ICR003007, CRE806, CRE1173, CRE346, 2580, CRK0001, and YDC667-1. Other isolation vere human illnesses. The clinical strains included AK-8, ICR003007, CRE806, CRE1173, sources for the strains in clade V included chicken (strains JS97 and L38), sprouts (strain CRE346, 2580, CRK0001, and YDC667-1. Other isolation sources for the strains in clade C5.1), and coastal water (strain CW_LB-887), and a significant correlation was identified V included chicken (strains JS97 and L38), sprouts (strain C5.1), and coastal water (strain for clade V strains with sampling sources from diverse geographical locations, excluding CW_LB-887), and a significant correlation was identified for clade V strains with sampling the United States (Fisher’s exact test p-value of < 0.02). sources from diverse geographical locations, excluding the United States (Fisher’s exact test p-value of <0.02). Microorganisms 2023, 11, 2114 9 of 17 Microorganisms 2023, 11, x FOR PEER REVIEW 9 of 17 Additionally, the functional COG categories associated with the C. werkmanii accessory genomeAdditionally, are presented theinfunctional Figure 2. COG Whencategories examining associated the entire with the C. werkmanii pangenome, the average accessoryofgenome percentage are presented functional COG wasin85%, Figure and2. subsequent When examining the entire analysis of thepangenome, gene clusters thefor the average percentage of functional COG was 85%, and subsequent core, shell, and cloud genomes revealed the average percentage of the functional COGs to analysis of the gene clusters26.5%, be 91.5%, for theand core, shell, respectively. 57.8%, and cloud genomes revealedmetabolism Carbohydrate the average and percentage of the transport, cell wall functional COGs to be 91.5%, 26.5%, and 57.8%, respectively. Carbohydrate metabolism structure, translation, transcription, signal transduction, and general metabolism (amino and transport, cell wall structure, translation, transcription, signal transduction, and acids, lipids, and nucleotides) were detected as the principal orthologous groups of proteins general metabolism (amino acids, lipids, and nucleotides) were detected as the principal during the analysis orthologous groupsof of the coreduring proteins genome. These of the analysis orthologous the core genome.proteinThese groups in the core orthologous genome were associated with C. werkmanii traits that are needed protein groups in the core genome were associated with C. werkmanii traits that are needed for growth and survival, while for growth and survival, while the shell genome was found to be mainly composed of in the the shell genome was found to be mainly composed of genes that are implicated secretion system, genes that pilus assembly are implicated in theprotein, secretionand phage system, structural pilus assembly genes. protein,Whenand compared phage to the structural core genome, genes. When the compared analysis of thetocloud the coreandgenome, the analysis shell genomes of thesignificant showed cloud and shell variability. genomes showed significant By conducting a subsequent variability. analysis with the Anvi’o software version 7, the presence By conducting a subsequent of putative genes for strains belonging analysis to with eachthe Anvi’o clade wassoftware version identified and7,theythe presence were assigned of putative genes for strains belonging to each clade was to functional categories for the shell (Figure 3a) and cloud (Figure 3b) pangenomes. identified and they were The assigned to functional categories for the shell (Figure 3a) and cloud (Figure 3b) cellular process and signaling (brown-yellow) category, including genes with a proposed pangenomes. The cellular process and signaling (brown-yellow) category, including function in kinase signal transduction, accounted for about 26% to 39% of the genes in the genes with a proposed function in kinase signal transduction, accounted for about 26% to shell 39%pangenome. of the genesApproximately 25% to 32% in the shell pangenome. of the genes Approximately 25% were assigned to 32% of the to thiswere genes functional category assignedin the cloud to this pangenome. functional category Moreover, in the cloud thepangenome. informationMoreover, storage and processing (blue) the information category, storageincluding genes(blue) and processing with category, a proposed function including genesin transcription, with a proposed showed function theinlowest percentage, ranging transcription, showedfrom the14% lowestto percentage, 18% of theranginggenes in fromthe14%shell to genome 18% of theand 12% genes in to the21% of the shell genesgenome in theand 12% pangenome. cloud to 21% of the genes in the cloudthe Interestingly, pangenome. Interestingly, most variability wasthe most observed for variability metabolic was observed functions for metabolic (pink-red), such asfunctions (pink-red), the transport such as of lipids, the transport of amino carbohydrates, lipids, acids, andcarbohydrates, nucleotides. amino acids, and The results nucleotides. for this category The hadresults for this category the highest percentageshad the ofhighest 34 to 38% of the genes, as detected for the human strains in clade I and the environmental clade the percentages of 34 to 38% of the genes, as detected for the human strains in clade I and IV for the environmental clade IV for the cloud and shell genomes, respectively. In addition, a high cloud and shell genomes, respectively. In addition, a high percentage of uncharacterized percentage of uncharacterized genes (grey-black), ranging from 23–33%, were genes (grey-black), ranging from 23–33%, were predominantly identified in the cloud predominantly identified in the cloud genomes for strains belonging to all clades, and the genomes for strains cloud genomes forbelonging the human to all clades, strains in clade and the clade II and cloudVgenomes for thebyhuman were identified havingstrains in clade II and clade V were identified by having prophages, prophages, transposons, and over 20% of the genes associated with the mobilome, transposons, and over 20% of the corresponding genes associated with the mobilome, corresponding to those to those genetic elements that can confer movement within and among thegenetic elements that can confer movement different bacterialwithin genomes and[54]. among the different bacterial genomes [54]. (a) (b) Figure 3. Functional COG categories derived from the Citrobacter werkmanii pangenome analysis. Over thirty C. werkmanii genomes (Table 2) were analyzed using Anvi’o software version 7, and putative genes were assigned to functional clusters of the orthologous gene categories for each clade, as described in Figure 2, for the shell genome (a) and cloud genome (b). Figure 3. Functional COG categories derived from the Citrobacter werkmanii pangenome analysis. Over thirty C. werkmanii genomes (Table 2) were analyzed using Anvi’o software version 7, and putative genes were assigned to functional clusters of the orthologous gene categories for each clade, as described in Figure 2, for the shell genome (a) and cloud genome (b). Microorganisms 2023, 11, 2114 10 of 17 3.3. Virulence and Antimicrobial Resistance Gene Identification To gain insights about specific virulence and resistance markers in C. werkmanii, the 3.3. Virulence and Antimicrobial Resistance Gene Identification available genome data from a collection of strains recovered from various geographical To gain insights about specific virulence and resistance markers in C. werkmanii, the locations availableand sources genome datawere from aanalyzed collection (Table 2).recovered of strains Additional fromCitrobacter species that are various geographical implicated locations in andhuman sourcesillness, the C. koseri were analyzed strain (Table ATCC BAA-895 2). Additional andspecies Citrobacter the C. that freundii are strain FDAARGOS_549, implicated in human illness, the C. koseri strain ATCC BAA-895 and the C. freundii strainwith a were added to the comparative analysis. A heatmap phylogenetic FDAARGOS_549,tree was wereconstructed added to theby indicatinganalysis. comparative the presence or absence A heatmap of targeted genes with a phylogenetic tree that arewas constructed implicated by indicating in virulence andthe presence or resistance antimicrobial absence of (Figure targeted4). genes that are implicated in virulence and antimicrobial resistance (Figure 4). Figure 4. Virulence Figure 4. Virulenceand andantimicrobial resistance antimicrobial resistance determinants determinants inexamined in the the examined Citrobacter Citrobacter werkmanii werkmanii strains. AA strains. heatmap, heatmap,indicating the presence indicating the presenceoror absence absence of virulence of virulence (top (top panel)panel) and antimicrobial and antimicrobial resistance genes (bottom panel), was constructed, coupled with a maximum-likelihood resistance genes (bottom panel), was constructed, coupled with a maximum-likelihood phylogenetic phylogenetic tree tree for the C. for the strains werkmanii C. werkmanii strains belonging to thebelonging to the various clades, various which clades, are shown which are in distinct shown in colors, distinct colors, as described in Figure 2, along with the clinical strains of C. koseri ATCC as described in Figure 2, along with the clinical strains of C. koseri ATCC BAA-895 and C. freundii BAA-895 andFDAARGOS_549. C. freundii FDAARGOS_549. Gene presence in the heatmap is indicated by the Gene presence in the heatmap is indicated by the red color and gene absence isred color and gene absence indicated byisthe indicated by the white color. whitevalues Bootstrap color.greater Bootstrap than values 85% aregreater thanlight shown with 85%blue arecircles shown with light over blue circles over the phylogenetic the phylogenetic tree nodes. tree nodes. Microorganisms 2023, 11, 2114 11 of 17 The analysis revealed that virulence traits related to the flagellum, type I, and type Va secretion systems were detected consistently in all of the examined Citrobacter strains. Interestingly, the phylogenetic analysis revealed the presence of mandatory secretion system-specific genes implicated in the type II and type Vb secretion system for the strains in clade V (Figure 4, top panel), comprising the largest number of examined C. werkmanii strains that were predominantly recovered from the clinical samples. The virulence map- ping analysis also revealed that more than 55% of the examined C. werkmanii strains were observed to harbor genes related to the type VI-class 1 secretion system. Additionally, dis- tinct profiles were also observed when examining the presence of the adhesin curli fimbrae csgF gene, the entS enterobactin gene, and the siderophore iroBCDN operon among the strains belonging to the different clades (Figure 4, top panel). In contrast, the adhesin genes csgB and csgE and the iron uptake genes entB, entE, and fepC were present in all examined C. werkmanii strains. By using the VirulenceFinder computational tool, genes implicated in the type VI secretion system were identified for C. werkmanii strains belonging to all five clades. Although a larger number of virulence genes were detected for the clinically relevant species, C. koseri, a total of 16 virulence genes were commonly detected for the C. werkmanii strains recovered from clinical, food, and environmental sources at various geographical locations. Multidrug resistance was observed in the examined C. werkmanii strains, and resistance to beta-lactamase, efflux pumps, peptide antibiotics, fluoroquinolone, and aminocoumarin were commonly observed in this study (Figure 4, bottom panel). With the exception of the strains FDAARGOS_616 and CRK0001, the remaining C. werkmanii strains that were clustering in clade V, which were predominantly from clinical sources, were found to have a larger number of genes associated with antimicrobial resistance in various drug classes when compared to the other examined C. werkmanii strains, as well as the clinical C. freundii FDAARGOS_549 and C. koseri ATCC BAA-895 reference strains (Figure 4, bottom panel). In particular, distinct antimicrobial gene profiles were detected among the C. werkmanii strains belonging to clade V for aminoglycoside {aac(3)-lle, ant(200 )-la, ant(300 )-lla, and ant(300 )-la}, beta-lactamases (blaOXA-10 , blaOXA-2 , blaOXA-370, and blaOXA-48 ), fluoroquinolone (qnrA1, qnrB11, qnrB12, qnrB17, and qnrB20), phenicol (catl, catll, cmlA1, and cmlA5), trimethoprim (dfrA1, dfrA12, dfrA14, dfrA17, dfrA19, dfrA23, and dfrA27), macrolide (mphA and mphE), sulfonamide (sul1, sul2, and sul3) and tetracycline (tetA, tetB, tetC, tetD, and tetM). The current analysis revealed that C. werkmanii strains belonging to clade V that were recovered from clinical sources showed the most extensive antimicrobial resistance profiles. In particular, a significantly higher number of resistance genes was observed for C. werkmanii, which had a total of 27 genes, when compared to the clinically relevant C. koseri strain ATCC BAA-895, with a total of 21 genes. In summary, these findings indicated that the C. werkmanii strains harbor virulence determinant genes, conferring an ability regarding host adaptation and colonization as well as multidrug resistance to various high-priority agents, indicating significant potential for C. werkmanii as an emerging and opportunistic pathogen in humans. 4. Discussion In the present study, a size-exclusion ultrafiltration method, followed by genome se- quencing and phylogenomic analyses, revealed the isolation of novel strains of C. werkmanii, the strains LANIIA-031 and LANIIA-032, which were recovered from agricultural river water in an important region for the production of fresh produce in northwestern Mexico. The examined C. werkmanii strains, LANIIA-031 and LANIIA-032, were isolated from major agricultural river samples as part of a field survey study aimed at the isolation of Salmonella enterica [21]. Preliminary typing and biochemical standardized assays falsely identified the LANIIA-031 and LANIIA-032 strains as Salmonella. However, the use of whole-genome sequencing led to the accurate species identification of the LANIIA-031 and LANIIA-032 strains as Citrobacter. In previous studies, phylogenomic analyses have indi- cated that Salmonella and Citrobacter have similar environmental niches and a highly related Microorganisms 2023, 11, 2114 12 of 17 evolutionary history [55]. The genomic similarity between Salmonella and Citrobacter poses a challenge in food safety since current traditional methods may lead to the inaccurate speciation of these foodborne pathogens. These findings indicate that the use of high- resolution sequencing-based methods focused on genomic characterization enables the accurate species identification of the recovered bacterial pathogens from diverse samples relevant to food safety, environmental monitoring, and diagnostics testing. Pangenome analyses of a large collection of sequenced C. werkmanii strains from diverse sources and locations resulted in the grouping of the strains into five distinct clades. Based on the identification of a gene’s presence or absence, as well as single nucleotide polymorphisms, most of the strains clustered in clade IV were from environmental sources and those in clade V were predominantly from clinical sources. The pangenome analysis from this study showed that the examined C. werkmanii strains had genomic features and a similar number of genes that were associated with subsystem categories, a finding that was in agreement with previous reports characterizing the genomic structure for another environmental C. werkmanii isolate recovered from industrial spoilage samples [56]. The current analysis also revealed that C. werkmanii exhibited an open pangenome, and this finding was in agreement with the characteristics of other Enterobacteriaceae family members such as S. enterica and Escherichia coli [57,58]. The supporting evidence indicating an open pangenome for C. werkmanii would imply a genome in a species that is constantly exchanging genetic material and implies an ability to colonize multiple environments, as observed for other pathogenic and environmental bacterial species [59]. In agreement with previous studies examining the pangenome of other Citrobacter species [60], the core pangenome in C. werkmanii comprised the essential gene components associated with general metabolism, replication, cell wall structure, and transcription categories in C. werkmanii strains belonging to all clades for this opportunistic pathogen. Although small variations between the functional gene classification among the strains were observed, the functional gene category analysis revealed that the shell and cloud genomes showed the highest variability in terms of the number of genes for cellular metabolism, including carbohydrate lipid, nucleotide, and amino acid transport, as well as secondary metabolite biosynthesis. In particular, carbohydrate metabolism genes were detected as the main accessory features, resulting in additional gene copies in the shell and cloud genomes for C. werkmanii strains belonging to all clades. In particular, genes with an important role in galactose metabolism were identified that were coding for galactose mutarotase and beta- galactosidase, which are involved in substrate conversion and degradation [61]. Previous studies have shown that galactose is an important carbohydrate for the lipopolysaccharide composition of Citrobacter species [62,63]. Interestingly, the cloud genomes for the human strains in clades II and V had a large number of prophages, transposons, and mobilome-associated genes. The acquisition of mobile elements is considered an important evolutionary mechanism for bacterial pathogens to adapt and colonize many reservoirs, as well as to enhance the virulence potential [64–66]. Additionally, the presence of phage protein tail and capsid structural genes were found as part of both the shell and cloud genomes. Phages have been described as a common bacterial evolution mechanism involving horizontal gene transfer [54,66]. Moreover, the findings from this study indicated a high number of genes of unknown function or uncharacterized genes, which were also in the cloud and shell genomes. These results indicate the current challenges in the annotation of microbial genomes, supporting the need for improved computational pipelines to integrate uncharacterized genes in an evolutionary, biotechnological, and ecological framework as the volumes of genomic data continue to be deposited in public repositories for the improved characterization of a particular bacterial species [67,68]. An analysis of virulence determinants showed that C. werkmanii strains in this study harbor genes that are implicated in secretion, as with the flagellum, type I, and type Va secretion systems. To our knowledge, the present study is the first one to identify specific clades associated with genetic traits linked to the secretion system in this species. Microorganisms 2023, 11, 2114 13 of 17 Moreover, genes related to the type II and type Vb secretion systems were found among the C. werkmanii strains that belonged to clade V, which were recovered predominantly from clinical sources. The type II secretion system is implicated in protein effector translocation across the outer membrane in pathogenic bacteria [69], while the type Vb secretion has the characteristic structure of a two-partner secretion pathway, with a role in the adhesion to receptors in mammalian cells [70]. Moreover, the identified type II secretion system genes in this study have a known role in the translocation of virulence factors, toxins, and enzymes across the cell’s outer membrane, as described for other human pathogens such as E. coli, Vibrio cholerae, Chlamydia trachomatis, and Acinetobacter baumannii [71]. Finally, genes of the type VI secretion system were identified in the C. werkmanii strains examined in this study. In a previous report, the type VI secretion system in C. freundii has been shown to contribute to adhesion and cytotoxicity in host cells [72]. Interestingly, the expression of the type VI secretion system by bacterial pathogens has been considered to be an antagonist mechanism in the colonization of the human gut by the resident microbiota [73], and the presence of the type VI secretion system could potentially provide C. werkmanii with a selective competitive advantage during the colonization of the human host. Another important observation of the virulence categorization was that the curli fim- briae (Csg) and outer membrane protein A (OmpA) were detected in most of the examined C. werkmanii strains. In particular, curli has been well characterized in pathogenic E. coli and S. enterica and has a role in biofilm formation and in attachment to the extracellular matrix and colonization of the mammalian host cells [74]. Moreover, recent studies on the C. werkmanii strain BF-6 have demonstrated that OmpA plays an important role in the regulation of multiple virulence phenotypes, including biofilm formation on distinct surfaces, swimming motility, metal ion responses, and resistance to biocides [75]. Finally, the virulence typing analyses of C. werkmanii indicated the presence of the siderophore enterobactin genes for the acquisition of iron, an essential nutrient for bacteria growth and the increased virulence potential of Citrobacter [76,77]. Multidrug resistance to various antimicrobial agents was commonly observed among the examined C. werkmanii strains in this study. Several studies have proposed beta-lactam resistance as a common mechanism among Citrobacter species [78]. Resistance to the fluoro- quinolone and aminocoumarin agents has previously been detected for other species of Citrobacter [60,77,79,80]. Interestingly, the present study demonstrated that the C. werkmanii strains that belong to clade V, which were predominantly from clinical sources, showed the most extensive antimicrobial resistance profiles. These observed multidrug resistances included various classes of antimicrobials, such as aminoglycoside, beta-lactamase, flu- oroquinolone, phenicol, trimethoprim, macrolide, sulfonamide, and tetracycline. Due to the mortality rate associated with ineffective therapy in the treatment of Citrobacter infections [13], these observations highlight the importance of the resistome characteriza- tion of this opportunistic bacterial species. The global dissemination of multidrug-resistant enteric pathogens, including Citrobacter, has resulted in a high-priority action item under the One Health approach by integrating health solutions in various interfaces related to hu- mans, animals, and the environment [81,82]. The spread of multidrug-resistant nosocomial infections to other niches is a critical issue that is in need of the improved identification and surveillance of bacterial virulence mechanisms and agent resistance [81]. In summary, this study demonstrated that the use of next-generation sequencing platforms in conjunc- tion with bioinformatics allowed the identification in C. werkmanii of markers that are linked to virulence and antimicrobial resistance, indicating the potential of this species as an emerging and opportunistic pathogen in humans. Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/microorganisms11082114/s1, Table S1: List of Citrobacter reference strains used in the present study; Table S2: Percentage of the average nucleotide identity for the Citrobacter strains examined in the present study. Microorganisms 2023, 11, 2114 14 of 17 Author Contributions: Conceptualization, J.R.A.-S., I.F.V.-L. and C.C.; methodology, J.R.A.-S., B.Q., R.P.-A. and B.G.L.; software, J.R.A.-S., J.A.O.-M., J.M.-U. and B.G.L.; validation, J.R.A.-S., I.F.V.-L. and B.L.; formal analysis, J.R.A.-S., B.Q., J.A.O.-M. and B.G.L.; investigation, J.R.A.-S., J.A.O.-M., R.P.-A. and B.G.L.; resources, J.R.A.-S., B.Q., J.M.-U. and C.C.; data curation, J.R.A.-S. and B.G.L.; writing—original draft preparation, J.R.A.-S., B.Q. and B.G.L.; writing—review and editing, J.R.A.-S., B.Q., B.G.L. and C.C.; visualization, J.R.A.-S., J.M.-U. and C.C.; supervision, J.R.A.-S., B.Q. and C.C.; project administration, C.C.; funding acquisition, B.Q. and C.C. All authors have read and agreed to the published version of the manuscript. Funding: This material was based in part upon work supported by the Laboratorio Nacional para la Investigación en Inocuidad Alimentaria (LANIIA) at Centro de Investigación y Desarrollo A. C. (CIAD) in Culiacán, Sinaloa and by the United States Department of Agriculture (USDA), Agricultural Research Service (ARS), CRIS Project Number 2030-42000-055-00D. Data Availability Statement: The complete genome sequences of the C. werkmanii strains of LANIIA- 031 and LANIIA-032 can be found under the GenBank accession numbers JAUJUK000000000 and JAUJUL000000000, respectively, with the BioProject accession number PRJNA992637 in the National Center for Biotechnology Information (NCBI) BioProject database (https://www.ncbi.nlm.nih.gov/ bioproject/ accessed on 15 July 2023). Acknowledgments: The authors would like to thank Irvin González López (CIAD, Culiacán, Sinaloa, México) for coordinating the survey of agricultural river water in the Culiacan Valley, Sinaloa, Mexico, and Maria Brandl (USDA-ARS, Albany, CA, USA) for providing critical comments on the manuscript. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. USDA is an equal opportunity provider and employer. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. References 1. Rogers, L.; Power, K.; Gaora, P.Ó.; Fanning, S. Escherichia coli and other Enterobacteriaceae: Occurrence and detection. In Encyclopedia of Food and Health; Caballero, B., Finglas, P.M., Toldrá, F., Eds.; Academic Press: Oxford, UK, 2016; pp. 545–551. 2. Kus, J.V. Infections due to Citrobacter and Enterobacter. In Reference Module in Biomedical Sciences; Elsevier: Amsterdam, The Netherlands, 2014. 3. Janda, J.M.; Abbott Sharon, L. The changing face of the family Enterobacteriaceae (Order: “Enterobacterales”): New members, taxonomic issues, geographic expansion, and new diseases and disease syndromes. Clin. Microbiol. Rev. 2021, 34, e00174-20. [CrossRef] [PubMed] 4. Jung, W.J.; Kim, H.J.; Giri, S.S.; Kim, S.G.; Kim, S.W.; Kang, J.W.; Kwon, J.; Lee, S.B.; Oh, W.T.; Jun, J.W.; et al. Citrobacter tructae sp. nov. isolated from kidney of diseased rainbow trout (Oncorhynchus mykiss). Microorganisms 2021, 9, 275. [CrossRef] [PubMed] 5. Oberhettinger, P.; Schüle, L.; Marschal, M.; Bezdan, D.; Ossowski, S.; Dörfel, D.; Vogel, W.; Rossen, J.W.; Willmann, M.; Peter, S. Description of Citrobacter cronae sp. nov., isolated from human rectal swabs and stool samples. Int. J. Syst. Evol. Microbiol. 2020, 70, 2998–3003. [CrossRef] 6. Ribeiro, T.G.; Izdebski, R.; Urbanowicz, P.; Carmeli, Y.; Gniadkowski, M.; Peixe, L. Citrobacter telavivum sp. nov. with chromosomal mcr-9 from hospitalized patients. Eur. J. Clin. Microbiol. Infect Dis. 2021, 40, 123–131. [CrossRef] 7. Ong, C.C.H.; Farhanah, S.; Linn, K.Z.; Tang, Y.W.; Poon, C.Y.; Lim, A.Y.; Tan, H.R.; Binte Hamed, N.H.; Huan, X.; Puah, S.H.; et al. Nosocomial infections among COVID-19 patients: An analysis of intensive care unit surveillance data. Antimicrob. Resist. Infect Control 2021, 10, 119. [CrossRef] 8. Pletz, M.W.; Wollny, A.; Dobermann, U.-H.; Rödel, J.; Neubauer, S.; Stein, C.; Brandt, C.; Hartung, A.; Mellmann, A.; Trommer, S.; et al. A nosocomial foodborne outbreak of a VIM carbapenemase-expressing Citrobacter freundii. Clin. Infect Dis. 2018, 67, 58–64. [CrossRef] [PubMed] 9. Ranjan, K.P.; Ranjan, N. Citrobacter: An emerging health care associated urinary pathogen. Urol. Ann. 2013, 5, 313. [CrossRef] 10. Tchidjou, H.K.; Romeo, B. Infant case of co-infection with SARS-CoV-2 and Citrobacter koseri urinary infection. J. Trop. Pediatr. 2021, 67, fmaa032. [CrossRef] 11. Anderson, M.T.; Mitchell, L.A.; Zhao, L.; Mobley, H.L.T. Citrobacter freundii fitness during bloodstream infection. Sci. Rep. 2018, 8, 11792. [CrossRef] [PubMed] 12. Ariza-Prota, M.A.; Pando-Sandoval, A.; García-Clemente, M.; Fernández, R.; Casan, P. Community-acquired pneumonia and empyema caused by Citrobacter koseri in an immunocompetent patient. Case Rep. Pulmonol. 2015, 2015, 670373. [CrossRef] 13. Pepperell, C.; Kus, J.V.; Gardam, M.A.; Humar, A.; Burrows, L.L. Low-virulence Citrobacter species encode resistance to multiple antimicrobials. Antimicrob. Agents Chemother. 2002, 46, 3555–3560. [CrossRef] [PubMed] Microorganisms 2023, 11, 2114 15 of 17 14. Gardam, M.A.; Burrows, L.L.; Kus, J.V.; Brunton, J.; Low, D.E.; Conly, J.M.; Humar, A. Is surveillance for multidrug-resistant Enterobacteriaceae an effective infection control strategy in the absence of an outbreak? J. Infect. Dis. 2002, 186, 1754–1760. [CrossRef] [PubMed] 15. Hammerum, A.M.; Hansen, F.; Nielsen, H.L.; Jakobsen, L.; Stegger, M.; Andersen, P.S.; Jensen, P.; Nielsen, T.K.; Hansen, L.H.; Hasman, H.; et al. Use of WGS data for investigation of a long-term NDM-1-producing Citrobacter freundii outbreak and secondary in vivo spread of blaNDM-1 to Escherichia coli, Klebsiella pneumoniae and Klebsiella oxytoca. J. Antimicrob. Chemother. 2016, 71, 3117–3124. [CrossRef] 16. Jones, M.E.; Karlowsky, J.A.; Draghi, D.C.; Thornsberry, C.; Sahm, D.F.; Nathwani, D. Epidemiology and antibiotic susceptibility of bacteria causing skin and soft tissue infections in the USA and Europe: A guide to appropriate antimicrobial therapy. Int. J. Antimicrob. Agents 2003, 22, 406–419. [CrossRef] [PubMed] 17. Lee, R.; Choi, S.M.; Jo, S.J.; Lee, J.; Cho, S.Y.; Kim, S.H.; Lee, D.G.; Jeong, H.S. Clinical characteristics and antimicrobial susceptibility trends in Citrobacter bacteremia: An 11-year single-center experience. Infect. Chemother. 2019, 51, 1–9. [CrossRef] 18. Mohanty, S.; Singhal, R.; Sood, S.; Dhawan, B.; Kapil, A.; Das, B.K. Citrobacter infections in a tertiary care hospital in Northern India. J. Infect. 2007, 54, 58–64. [CrossRef] 19. Nayar, R.; Shukla, I.; Sultan, A. Epidemiology, prevalence and identification of Citrobacter species in clinical specimens in a tertiary care hospital in India. Int. J. Sci. Res. Publ. 2014, 4, 1–6. 20. Parvez, S.; Khan, A.U.; Kaur, G.; Barakat, M.; Ortet, P.; Mayilraj, S. An insight into the genome of extensively drug-resistant and uropathogenic Citrobacter werkmanii. J. Glob. Antimicrob. Resist. 2020, 22, 785–791. [CrossRef] 21. González-López, I.; Medrano-Félix, J.A.; Castro-Del Campo, N.; López-Cuevas, O.; González-Gómez, J.P.; Valdez-Torres, J.B.; Aguirre-Sánchez, J.R.; Martínez-Urtaza, J.; Gómez-Gil, B.; Lee, B.G.; et al. Prevalence and genomic diversity of Salmonella enterica recovered from river water in a major sgricultural region in Northwestern Mexico. Microorganisms 2022, 10, 1214. [CrossRef] [PubMed] 22. Aguirre-Sánchez, J.R.; Valdez-Torres, J.B.; del Campo, N.C.; Martínez-Urtaza, J.; del Campo, N.C.; Lee, B.G.; Quiñones, B.; Chaidez-Quiroz, C. Phylogenetic group and virulence profile classification in Escherichia coli from distinct isolation sources in Mexico. Infect. Genet. Evol. 2022, 106, 105380. [CrossRef] 23. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. Available online: http://www.bioinformatics. babraham.ac.uk/projects/fastqc (accessed on 1 May 2021). 24. Krueger, F.; Trim Galore. A Wrapper Tool around Cutadapt FastQC to Consistently Apply QUAL Adapt Trimming to FastQ Files. Available online: https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ (accessed on 1 May 2021). 25. Bushnell, B. BBTools Software Package. Available online: Sourceforge.net/projects/bbmap/ (accessed on 1 May 2021). 26. Coil, D.; Jospin, G.; Darling, A.E. A5-miseq: An updated pipeline to assemble microbial genomes from Illumina MiSeq data. Bioinformatics 2015, 31, 587–589. [CrossRef] 27. Seemann, T. Prokka: Rapid prokaryotic genome annotation. Bioinformatics 2014, 30, 2068–2069. [CrossRef] 28. Brettin, T.; Davis, J.J.; Disz, T.; Edwards, R.A.; Gerdes, S.; Olsen, G.J.; Olson, R.; Overbeek, R.; Parrello, B.; Pusch, G.D.; et al. RASTtk: A modular and extensible implementation of the RAST algorithm for building custom annotation pipelines and annotating batches of genomes. Sci. Rep. 2015, 5, 8365. [CrossRef] [PubMed] 29. Jia, B.; Raphenya, A.R.; Alcock, B.; Waglechner, N.; Guo, P.; Tsang, K.K.; Lago, B.A.; Dave, B.M.; Pereira, S.; Sharma, A.N.; et al. CARD 2017: Expansion and model-centric curation of the comprehensive antibiotic resistance database. Nucleic Acids Res. 2017, 45, D566–D573. [CrossRef] [PubMed] 30. Liu, B.; Zheng, D.; Jin, Q.; Chen, L.; Yang, J. VFDB 2019: A comparative pathogenomic platform with an interactive web interface. Nucleic Acids Res. 2019, 47, D687–D692. [CrossRef] 31. Seemann, T. ABRicate. Available online: https://github.com/tseemann/abricate (accessed on 1 May 2021). 32. Letunic, I.; Bork, P. Interactive Tree Of Life (iTOL) v4: Recent updates and new developments. Nucleic Acids Res. 2019, 47, W256–W259. [CrossRef] 33. Abby, S.S.; Néron, B.; Ménager, H.; Touchon, M.; Rocha, E.P.C. MacSyFinder: A program to mine genomes for molecular dystems with an spplication to CRISPR-Cas systems. PLoS ONE 2014, 9, e110726. [CrossRef] 34. Abby, S.S.; Rocha, E.P.C. Identification of protein secretion systems in bacterial genomes using MacSyFinder. In Bacterial Protein Secretion Systems: Methods and Protocols; Journet, L., Cascales, E., Eds.; Springer New York: New York, NY, USA, 2017; pp. 1–21. 35. Clermont, D.; Motreff, L.; Passet, V.; Fernandez, J.-C.; Bizet, C.; Brisse, S. Multilocus sequence analysis of the genus Citrobacter and description of Citrobacter pasteurii sp. nov. Int. J. Syst. Evol. Microbiol. 2015, 65, 1486–1490. [CrossRef] [PubMed] 36. Das, S.; Dash, H.R.; Mangwani, N.; Chakraborty, J.; Kumari, S. Understanding molecular identification and polyphasic taxonomic approaches for genetic relatedness and phylogenetic relationships of microorganisms. J. Microbiol. Methods 2014, 103, 80–100. [CrossRef] [PubMed] 37. Kumar, S.; Stecher, G.; Li, M.; Knyaz, C.; Tamura, K. MEGA X: Molecular Evolutionary Genetics Analysis across Computing Platforms. Mol. Biol. Evol. 2018, 35, 1547–1549. [CrossRef] [PubMed] 38. Saitou, N.; Nei, M. The neighbor-joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 1987, 4, 406–425. [CrossRef] [PubMed] 39. Tamura, K.; Nei, M.; Kumar, S. Prospects for inferring very large phylogenies by using the neighbor-joining method. Proc. Natl. Acad. Sci. USA 2004, 101, 11030–11035. [CrossRef] Microorganisms 2023, 11, 2114 16 of 17 40. Chun, J.; Oren, A.; Ventosa, A.; Christensen, H.; Arahal, D.R.; da Costa, M.S.; Rooney, A.P.; Yi, H.; Xu, X.W.; De Meyer, S.; et al. Proposed minimal standards for the use of genome data for the taxonomy of prokaryotes. Int. J. Syst. Evol. Microbiol. 2018, 68, 461–466. [CrossRef] [PubMed] 41. Meier-Kolthoff, J.P.; Carbasse, J.S.; Peinado-Olarte, R.L.; Göker, M. TYGS and LPSN: A database tandem for fast and reliable genome-based classification and nomenclature of prokaryotes. Nucleic Acids Res. 2022, 50, D801–D807. [CrossRef] [PubMed] 42. Richter, M.; Rossello-Mora, R.; Oliver Glockner, F.; Peplies, J. JSpeciesWS: A web server for prokaryotic species circumscription based on pairwise genome comparison. Bioinformatics 2016, 32, 929–931. [CrossRef] [PubMed] 43. Bertelli, C.; Laird, M.R.; Williams, K.P.; Simon Fraser University Research Computing Group; Lau, B.Y.; Hoad, G.; Winsor, G.L.; Brinkman, F.S.L. IslandViewer 4: Expanded prediction of genomic islands for larger-scale datasets. Nucleic Acids Res. 2017, 45, W30–W35. [CrossRef] 44. Arndt, D.; Grant, J.R.; Marcu, A.; Sajed, T.; Pon, A.; Liang, Y.; Wishart, D.S. PHASTER: A better, faster version of the PHAST phage search tool. Nucleic Acids Res. 2016, 44, W16–W21. [CrossRef] 45. Darling, A.E.; Mau, B.; Perna, N.T. progressiveMauve: Multiple genome alignment with gene gain, loss and rearrangement. PLoS ONE 2010, 5, e11147. [CrossRef] 46. Page, A.J.; Cummins, C.A.; Hunt, M.; Wong, V.K.; Reuter, S.; Holden, M.T.G.; Fookes, M.; Falush, D.; Keane, J.A.; Parkhill, J. Roary: Rapid large-scale prokaryote pan genome analysis. Bioinformatics 2015, 31, 3691–3693. [CrossRef] 47. Eren, A.M.; Esen, Ö.C.; Quince, C.; Vineis, J.H.; Morrison, H.G.; Sogin, M.L.; Delmont, T.O. Anvi’o: An advanced analysis and visualization platform for ‘omics data. PeerJ 2015, 3, e1319. [CrossRef] 48. Treangen, T.J.; Ondov, B.D.; Koren, S.; Phillippy, A.M. The Harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes. Genome Biol. 2014, 15, 524. [CrossRef] 49. Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 2014, 30, 1312–1313. [CrossRef] [PubMed] 50. Hadfield, J.; Croucher, N.J.; Goater, R.J.; Abudahab, K.; Aanensen, D.M.; Harris, S.R. Phandango: An interactive viewer for bacterial population genomics. Bioinformatics 2018, 34, 292–293. [CrossRef] [PubMed] 51. Tettelin, H.; Masignani, V.; Cieslewicz, M.J.; Donati, C.; Medini, D.; Ward, N.L.; Angiuoli, S.V.; Crabtree, J.; Jones, A.L.; Durkin, A.S.; et al. Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: Implications for the microbial “pan-genome”. Proc. Natl. Acad. Sci. USA 2005, 102, 13950–13955. [CrossRef] [PubMed] 52. Mehta, C.R.; Patel, N.R. Algorithm 643. FEXACT: A FORTRAN subroutine for Fisher’s exact test on unordered r × c contingency tables. ACM Trans. Math Softw. 1986, 12, 154–161. [CrossRef] 53. Goris, J.; Konstantinidis, K.T.; Klappenbach, J.A.; Coenye, T.; Vandamme, P.; Tiedje, J.M. DNA-DNA hybridization values and their relationship to whole-genome sequence similarities. Int. J. Syst. Evol. Microbiol. 2007, 57, 81–91. [CrossRef] 54. Smets, B.F.; Barkay, T. Horizontal gene transfer: Perspectives at a crossroads of scientific disciplines. Nat. Rev. Microbiol. 2005, 3, 675–678. [CrossRef] [PubMed] 55. Pilar, A.V.C.; Petronella, N.; Dussault, F.M.; Verster, A.J.; Bekal, S.; Levesque, R.C.; Goodridge, L.; Tamber, S. Similar yet different: Phylogenomic analysis to delineate Salmonella and Citrobacter species boundaries. BMC Genom. 2020, 21, 377. [CrossRef] [PubMed] 56. Zhou, G.; Peng, H.; Wang, Y.-s.; Huang, X.-m.; Xie, X.-b.; Shi, Q.-s. Complete genome sequence of Citrobacter werkmanii strain BF-6 isolated from industrial putrefaction. BMC Genom. 2017, 18, 765. [CrossRef] 57. Tantoso, E.; Eisenhaber, B.; Kirsch, M.; Shitov, V.; Zhao, Z.; Eisenhaber, F. To kill or to be killed: Pangenome analysis of Escherichia coli strains reveals a tailocin specific for pandemic ST131. BMC Biol. 2022, 20, 146. [CrossRef] 58. Laing, C.R.; Whiteside, M.D.; Gannon, V.P.J. Pan-genome analyses of the species Salmonella enterica, and identification of genomic markers predictive for species, subspecies, and serovar. Front Microbiol. 2017, 8, 1345. [CrossRef] [PubMed] 59. Vernikos, G.; Medini, D.; Riley, D.R.; Tettelin, H. Ten years of pan-genome analyses. Curr. Opin. Microbiol. 2015, 23, 148–154. [CrossRef] [PubMed] 60. Yuan, C.; Yin, Z.; Wang, J.; Qian, C.; Wei, Y.; Zhang, S.; Jiang, L.; Liu, B. Comparative genomic analysis of Citrobacter and key genes essential for the pathogenicity of Citrobacter koseri. Front Microbiol. 2019, 10, 2774. [CrossRef] 61. Thoden, J.B.; Kim, J.; Raushel, F.M.; Holden, H.M. The catalytic mechanism of galactose mutarotase. Protein Sci. 2003, 12, 1051–1059. [CrossRef] 62. Kocharova, N.A.; Knirel, Y.A.; Stanislavsky, E.S.; Kholodkova, E.V.; Lugowski, C.; Jachymek, W.; Romanowska, E. Structural and serological studies of lipopolysaccharides of Citrobacter O35 and O38 antigenically related to Salmonella. FEMS Microbiol. Immunol. 1996, 13, 1–8. [CrossRef] 63. Raff, R.A.; Wheat, R.W. Carbohydrate composition of the phenol-soluble lipopolysaccharides of Citrobacter freundii. J. Bacteriol. 1968, 95, 2035–2043. [CrossRef] [PubMed] 64. Aguirre-Sánchez, J.R.; Ibarra-Rodriguez, J.R.; Vega-López, I.F.; Martínez-Urtaza, J.; Chaidez-Quiroz, C. Genomic signatures of adaptation to natural settings in non-typhoidal Salmonella enterica serovars Saintpaul, Thompson and Weltevreden. Infect. Genet. Evol. 2021, 90, 104771. [CrossRef] 65. Subedi, D.; Kohli, G.S.; Vijay, A.K.; Willcox, M.; Rice, S.A. Accessory genome of the multi-drug resistant ocular isolate of Pseudomonas aeruginosa PA34. PLoS ONE 2019, 14, e0215038. [CrossRef] 66. Colavecchio, A.; Cadieux, B.; Lo, A.; Goodridge, L.D. Bacteriophages contribute to the spread of antibiotic resistance genes among foodborne pathogens of the Enterobacteriaceae family—A review. Front Microbiol. 2017, 8, 1108. [CrossRef] Microorganisms 2023, 11, 2114 17 of 17 67. Miller, D.; Stern, A.; Burstein, D. Deciphering microbial gene function using natural language processing. Nat. Commun. 2022, 13, 5731. [CrossRef] 68. Vanni, C.; Schechter, M.S.; Acinas, S.G.; Barberán, A.; Buttigieg, P.L.; Casamayor, E.O.; Delmont, T.O.; Duarte, C.M.; Eren, A.M.; Finn, R.D.; et al. Unifying the known and unknown microbial coding sequence space. eLife 2022, 11, e67667. [CrossRef] [PubMed] 69. Korotkov, K.V.; Sandkvist, M.; Hol, W.G.J. The type II secretion system: Biogenesis, molecular architecture and mechanism. Nat. Rev. Microbiol. 2012, 10, 336–351. [CrossRef] [PubMed] 70. Meuskens, I.; Saragliadis, A.; Leo, J.C.; Linke, D. Type V secretion systems: An overview of passenger domain functions. Front Microbiol. 2019, 10, 1163. [CrossRef] [PubMed] 71. Chernyatina, A.A.; Low, H.H. Core architecture of a bacterial type II secretion system. Nat. Commun. 2019, 10, 5437. [CrossRef] 72. Liu, L.; Hao, S.; Lan, R.; Wang, G.; Xiao, D.; Sun, H.; Xu, J. The type VI secretion system modulates flagellar gene expres- sion and secretion in Citrobacter freundii and contributes to adhesion and cytotoxicity to host cells. Infect. Immun. 2015, 83, 2596–2604. [CrossRef] 73. Chen, C.; Yang, X.; Shen, X. Confirmed and potential roles of bacterial T6SSs in the intestinal ecosystem. Front. Microbiol. 2019, 10, 1484. [CrossRef] 74. Barnhart, M.M.; Chapman, M.R. Curli biogenesis and function. Annu. Rev. Microbiol. 2006, 60, 131–147. [CrossRef] 75. Zhou, G.; Wang, Y.S.; Peng, H.; Li, S.J.; Sun, T.L.; Shen, P.F.; Xie, X.B.; Shi, Q.S. Roles of ompA of Citrobacter werkmanii in bacterial growth, biocide resistance, biofilm formation and swimming motility. Appl. Microbiol. Biotechnol. 2021, 105, 2841–2854. [CrossRef] 76. Kortman, G.A.; Boleij, A.; Swinkels, D.W.; Tjalsma, H. Iron availability increases the pathogenic potential of Salmonella typhimurium and other enteric pathogens at the intestinal epithelial interface. PLoS ONE 2012, 7, e29968. [CrossRef] 77. Qin, J.; Zhao, Y.; Wang, A.; Chi, X.; Wen, P.; Li, S.; Wu, L.; Bi, S.; Xu, H. Comparative genomic characterization of multidrug- resistant Citrobacter spp. strains in Fennec fox imported to China. Gut Pathog. 2021, 13, 59. [CrossRef] 78. Liu, L.; Lan, R.; Liu, L.; Wang, Y.; Zhang, Y.; Wang, Y.; Xu, J. Antimicrobial resistance and cytotoxicity of Citrobacter spp. in Maanshan Anhui Province, China. Front Microbiol. 2017, 8, 1357. [CrossRef] 79. Aoyama, H.; Fujimaki, K.; Sato, K.; Fujii, T.; Inoue, M.; Hirai, K.; Mitsuhashi, S. Clinical isolate of Citrobacter freundii highly resistant to new quinolones. Antimicrob. Agents Chemother. 1988, 32, 922–924. [CrossRef] [PubMed] 80. Metri, B.C.; Jyothi, P.; Peerapur, B.V. Antibiotic resistance in Citrobacter spp. isolated from urinary tract infection. Urol. Ann. 2013, 5, 312. [PubMed] 81. Butler Mark, S.; Gigante, V.; Sati, H.; Paulin, S.; Al-Sulaiman, L.; Rex John, H.; Fernandes, P.; Arias Cesar, A.; Paul, M.; Thwaites Guy, E.; et al. Analysis of the clinical pipeline of treatments for drug-resistant bacterial infections: Despite progress, more action is needed. Antimicrob. Agents Chemother. 2022, 66, e01991-21. [CrossRef] 82. Velázquez-Meza, M.E.; Galarde-López, M.; Carrillo-Quiróz, B.; Alpuche-Aranda, C.M. Antimicrobial resistance: One Health approach. Vet World 2022, 15, 743–749. [CrossRef] [PubMed] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.