Proteus mirabilis causes complex urinary tract infections, which may lead to sepsis or systemic inflammatory response syndrome, with a fatality rate of 20%–50% ^[1]. P. mirabilis, the most common pathogen of the genus Proteus, is a member of the Enterobacteriaceae that causes urinary system infections, second only to Escherichia coli. P. mirabilis accounts for 3% of nosocomial infections ^[2]. With the widespread irregular use of antibiotics, P. mirabilis antibiotic resistance has increased. Multidrug-resistant P. mirabilis producing broad-spectrum β-lactamase was reported for the first time in Italy ^[3]. Subsequently, Poland, China, Japan, and other countries and regions have reported multidrug-resistant P. mirabilis ^[4-8], showing global spread. Due to the continuous increase in antibiotic resistance, the prevention and control of P. mirabilis infections is challenging. Therefore, high-resolution P. mirabilis typing is needed.

High-resolution sequence typing is important to assess the epidemiological links and identify possible sources of transmission during outbreaks. Currently, multilocus sequence typing (MLST) ^{[9, 10]} and pulsed-field gel electrophoresis (PFGE) ^[11-13] are the two most commonly used typing methods for outbreak investigations. However, their main disadvantages are the low resolution in MLST ^{[14, 15]} and difficulties in standardization of PFGE ^[16]. To the best of our knowledge, there are few studies on MLST of P. mirabilis, and most studies have focused on MLST of carbapenem-resistant Gram-negative clinical isolates containing P. mirabilis ^{[9, 10]}.

In recent years, with the continuous development of genomic epidemiology, the gene-by-gene (GbG) method has been advocated as an extension of MLST ^[17]. This method has the advantages of portability, scalability, and independence. It is regarded by PulseNet International as an important method for the identification of bacterial strains by high-throughput sequencing. The open source chewBBACA algorithm ^[18] provides a simple bioinformatics pipeline to construct cgMLST schemes of target strains. To type P. mirabilis globally and monitor the clonal groups (CGs) of potentially high-risk P. mirabilis, we used chewBBACA as a bioinformatics pipeline for GbG analysis to construct and verify the cgMLST scheme of P. mirabilis. Our results support the prevention and control of P. mirabilis infections.

Genomic sequences available as of September 10, 2021, for P. mirabilis (Enterobacteria) were downloaded from the National Center for Biotechnology Information (NCBI) genome sequence repository (http://www.ncbi.nlm.nih.gov/genome/). The downloaded sequences included 74 complete genomes and 679 whole genome shotgun sequences available as scaffolds (116), chromosomes (6), or contigs (557). After removing low-quality assemblies (gene assembly rate < 95% or contamination rate > 5%), all remaining genomes were included in the preliminary analysis.

In the present study, a chewBBACA suite was designed to create and evaluate the core genome GbG typing scheme, followed by allele calling in available P. mirabilis strains. ChewBBACA performs scheme creation and allele calls for complete or draft genomes generated by de novo assembly. Before starting, we created a training file in Prodigal v2.6.3 ^[19] using the reference genome HI4320 of P. mirabilis. The genome sequence of P. mirabilis HI4320 (RefSeq assembly accession: GCF_000069965.1) was used only as a reference genome to predict the whole genome MLST (wgMLST) loci and was then removed from further analysis.

A total of 72 complete genome sequences were annotated in the first step. In this step, the algorithm defined the coding sequences (CDSs) of each genome, and then all CDSs in the genome were compared in a pairwise manner to generate a single FASTA file containing the selected CDSs. Then, a two-step evaluation process was conducted to identify all CDSs in the genomes. First, CDSs that have the same sequence but are smaller than other CDSs were removed, and the larger CDSs were retained. Second, the remaining CDSs were gathered in unique loci by performing an all-against-all BLASTP search and calculating the blast score ratio (BSR) ^[20]. CDSs with a pairwise BSR of > 0.6 were considered to be alleles at the same locus, and the larger allele (in bp) was retained. The remaining wgMLST loci were used to define the cgMLST scheme. We selected candidate loci from 100% (threshold) of the available complete genomes (72 genomes) to define our cgMLST scheme. The second step was validation of candidate loci based on a total of 635 unfinished P. mirabilis genomes. We reduced stringency and maintained candidate loci shared by 99% of the isolates. This cutoff was chosen based on the quality of the tested genomes, since the sequences used in this step were incomplete genomes. ChewBBACA is available at https://github.com/B-UMMI/chewBBACA.

To define CGs in a non-arbitrary manner, we analyzed the distribution of the number of allelic mismatches (loci with different sequences) among the pairs of genomes. We performed a single-link algorithm designed to classify isolates from cgMLST data. Similar to the classification of classical MLST data, the single-link algorithm was also used for cgMLST data to accurately classify isolates ^[21].

A minimum spanning tree (MST) was constructed using the MSTree method in GrapeTree (version 2.1) for the allelic profiles obtained for each isolate by the cgMLST scheme.

Based on our cgMLST scheme, MSTs for different geographical locations were built to investigate possible evolutionary relationships. We focused on the distribution of CGs in different countries or regions and examined the shared and unique CGs between them.

We used the ABRicate pipeline (https://github.com/tseemann/abricate) to annotate the antibiotic resistance genes (ARGs) and virulence genes in the P. mirabilis genome using the CARD (https://card.mcmaster.ca/) and VFDB (http://www.mgc.ac.cn/VFs/main.htm) databases. The annotation parameters for ARGs and virulence genes were gene coverage > 60% and identity > 70%.

As of September 10, 2021, 753 genome sequences were available at the NCBI genome sequence repository (http://www.ncbi.nlm.nih.gov/genome/). Of these genomes, 45 were not included due to the low assembly quality (Supplementary Table S1 available in www.besjournal.com). The genome sequences of 12 P. mirabilis strains were excluded during the cgMLST scheme validation process (Supplementary Table S1). According to our cgMLST scheme, reference genome HI4320 was not included. Finally, 695 P. mirabilis strains were included in the analysis according to our cgMLST scheme. Information regarding all strains is provided in Supplementary Table S2 (available in www.besjournal.com). As shown in Figure 1, the majority of P. mirabilis isolates were collected from humans (82.4%), most samples were collected in the USA or China (75.3%), urine and sputum samples were the most common (40.6%), and most samples were collected in 2018 or 2019 (54.8%).

Figure 1.  Characterization of the P. mirabilis isolates included in this study. (A) Hosts of isolated P. mirabilis. “Others” represents host types with only one strain. (B) Geographic distribution of available genomes. “Others” represents countries with less than five strains. (C) Isolation source distribution of isolates. “Others” represents sources with less than eight strains. (D) Temporal distribution of P. mirabilis isolates. “Others” represents years when less than two strains were collected.

A flowchart of the cgMLST scheme of P. mirabilis (https://github.com/B-UMMI/chewBBACA) is provided in Figure 2. In total, 72 complete genome sequences were retrieved. A total of 7,305 target loci were generated and annotated using the wgMLST dataset. Using the RemoveGenes operation to remove possible paralogous loci from the wgMLST scheme, we discarded 85 loci. We filtered out an additional 4,921 loci using the genome quality test (threshold 55). A total of 2,299 candidate loci for cgMLST were retained, which were present in 72 complete genome sequences (Supplementary Figure S1 available in www.besjournal.com). During the creation of the scheme, no genome sequence was removed.

Figure 2.  Flowchart describing the development of the cgMLST scheme for P. mirabilis using chewBBACA (https://github.com/B-UMMI/chewBBACA).

We used the 635 unfinished P. mirabilis genomes for the validation of the cgMLST scheme. The TestGenomeQuality operation was used to select loci present in a predetermined percentage (99%) of the analyzed strains from the cgMLST scheme (Supplementary Figure S2 available in www.besjournal.com). During the validation step, 12 unfinished genomes were removed. In the end, 1,842 cgMLST loci were selected. Following this step, a cgMLST scheme consisting of 1,842 gene targets was defined, covering 49.91% of the 3,690 open reading frames predicted for the reference strain HI4320.

Figure S2.  Number of loci and number of genomes at every exclusion threshold level for the 635 validation strains’ genomes and 72 completed genomes that were used to construct the cgMLST scheme

The allelic difference threshold analysis showed a high number of genome pairs with < 291 mismatches or > 1,345 mismatches, and almost no genome pairs had 291–1,345 mismatches (Figure 3). Finally, we defined the P. mirabilis cgMLST CGs as a set of cgMLST profiles having ≤ 291 allelic mismatches with at least one other member of the group. The final 205 CGs were identified based on single-linkage clustering of the cgMLST allele profile.

Figure 3.  Distribution of the number of paired allele mismatches in P. mirabilis (i.e., the number of different loci for a given pair of strains). The allelic mismatch cutoff value of 291 proposed for the CG definition is shown in red.

From 695 P. mirabilis strains that conform to the cgMLST scheme, 205 CG types were identified (Supplementary Table S2). Owing to the limitation of the MST visualization, Figure 4 only exhibits the CGs (ranked by the number of strains included in the top 30) distributed in clusters. As shown in Figure 4, CGs exhibited good clustering. In other words, most P. mirabilis strains have been well typed according to our cgMLST scheme.

Figure 4.  Minimum spanning tree (MST) of CGs generated using the 695 P. mirabilis strains based on the cgMLST scheme. Due to the limitation of the visualization of the MST, only the CG types containing the top 30 strains are displayed. “Others” represents the CG types containing less than five strains.

Based on our cgMLST scheme, a minimal generative tree of 695 P. mirabilis strains from different countries or regions showed no close relationships at the geographic level (Figure 5). Overall, it could be seen that P. mirabilis presented a clear global spread, and the root of the MST originated from the USA. However, there are differences between shared CGs and unique CGs in different countries or regions. In total, 46 shared CGs were distributed in 26 countries or regions. CG3 was present in at least 13 countries or regions, CG20 was present in seven countries, and CG7 was present in five countries (Supplementary Figure S3, available in www.besjournal.com). Moreover, 159 unique CGs were distributed in 16 countries (Supplementary Table S3, available in www.besjournal.com). The USA had 105 unique CGs, China had 10 unique CGs, and Denmark had 13 unique CGs. Overall, there were regional differences in the distribution of CGs.

Figure 5.  Minimal spanning tree analysis of cgMLST profiles of 695 P. mirabilis strains isolated from different countries or regions. The genome sequences were downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/ genome/). Colors indicate the geographic locations of the isolates.

In total, 64 virulence gene types were found in 695 P. mirabilis strains. Supplementary Table S4 (available in www.besjournal.com) presents the distribution of the 64 virulence gene types in the 205 CGs. The carrier rate of virulence genes cheB, cheY, flgG, flgH, fliI, fliN, flip, gmhA/lpcA, kdsA, lpxC, lpxD, and luxS was 100%. The carrier rate of virulence genes cheW, flgC, flhC, flhD, fliA, fliG, fliM, fliZ, galU, htpB, ompA, and rfaD was > 95%. It is worth noting that CG20 carries nine unique virulence types, namely, papC, papD, papE, papF, papG, papH, papI, papJ, and papK. These nine unique types of virulence genes were present only in CG20 and were P fimbrial operons that cause severe urinary tract infections ^[22-24].

The subtypes of 156 ARGs were found in 695 P. mirabilis strains. Supplementary Table S5 (available in www.besjournal.com) presents the distribution of the subtypes of 156 ARGs among the 205 CGs. The ARG CRP was present in all CGs. The ARGs H-NS, acrB, cpxA, and tet(J) were present in more than 95% of CGs of P. mirabilis. The total number of ARG subtypes present in different CGs ranged from six to 95. CG3, CG20, CG29, CG32, CG6, CG49, CG41, CG80, CG12, CG7, CG4, CG1, CG61, CG68, CG48, CG28, CG107, and CG11 not only have a large number of shared ARG subtypes, but also unique ARG subtypes. CG3 has 77 shared ARG subtypes and 18 unique ARG subtypes. CG20 has 37 shared ARG subtypes and five unique ARG subtypes.

The spread of multidrug-resistant P. mirabilis ^[4-8] is an urgent public health problem. Proper typing of P. mirabilis is important. MLST ^{[9, 10]} and PFGE ^[11-13] are almost incapable of high-resolution typing of P. mirabilis. CgMLST is based on comparison of bacterial whole genome sequences, and has proven to be a promising means for disease outbreak investigation, source tracing, surveillance of bacterial pathogens, and contamination control ^[25-27]. cgMLST is based on comparison of a complete or draft de novo genome assembly with a scheme consisting of a series of loci and a collection of related allele sequences. This is a suitable and unbiased means to identify possible clusters from a sample of the entire species, and does not require any outbreak-specific reference ^[28]. To the best of our knowledge, no previous studies have proposed a cgMLST scheme for P. mirabilis typing.

GbG approaches are becoming increasingly popular in bacterial typing ^{[29, 30]}, bacterial genome epidemiology ^[31], and outbreak investigation ^{[31, 32]}. However, the lack of open source software for scheme definition and allele calling has limited the widespread use of the cgMLST approach. Interestingly, the chewBBACA suite can address this limitation. The chewBBACA suite was designed to help users create and evaluate novel whole genome or core genome GbG schemes and make allelic calls to bacterial strains of interest ^[17]. We used chewBBACA to call alleles of P. mirabilis and created and validated the first cgMLST scheme for P. mirabilis. In the end, the proposed cgMLST scheme for P. mirabilis typed 100% (707/707) of the genomes available for this research considering at least 95% of the cgMLST target genes present. We proposed the public cgMLST scheme for P. mirabilis based on 72 complete genome sequences and validated the scheme using 1,842 target loci in 635 unfinished genome sequences. The proposed 1842-locus cgMLST scheme delineated the 695 P. mirabilis strains into 205 distinct CGs. CG3 (18.3%) was the most dominant, followed by CG2 (4.7%), CG17 (3.0%), and CG32 (2.4%). Several CGs showed interesting regional distribution characteristics.

The cgMLST approach, which extends the MLST concept to the core genome, has proven to be a useful high-resolution typing scheme in other bacterial species ^{[29, 30, 33, 34]}. In contrast to the genus-wide Leptospira core genome MLST for strain taxonomy ^[30], we established a cgMLST scheme at the species level. Currently, the genus Proteus consists of five species (P. mirabilis, P. vulgaris, P. penneri, P. myxofaciens, and P. hauseri), with P. mirabilis being the main pathogenic species ^[35]. P. mirabilis is generally more susceptible to antimicrobials than P. penneri, P. hauseri, and P. vulgaris ^{[35, 36]}. Therefore, a cgMLST scheme at the species level is necessary. Unlike the cgMLST scheme of de Sales RO et al. for Pseudomonas aeruginosa ^[29], our cgMLST scheme allows for analysis of ARGs and virulence genes. The identification of CGs of high-risk P. mirabilis is important for public health.

Virulence genes papC, papD, papE, papF, papG, papH, papI, papJ, and papK were only in CG20. These nine virulence genes encode adhesive pili of the chaperone-usher family and they form the important P. mirabilis pmf fimbrial operon. These virulence genes in P. mirabilis may be the cause of severe urinary tract infections ^{[22-24, 37-39]}. In addition, CG20 has 42 shared ARG subtypes and five unique ARG subtypes. ARGs AAC(3)-Id ^{[40, 41]}, aadA7 ^[42,43], dfrA15 ^[44,45], and vanRA ^[46] have been found in multidrug-resistant P. mirabilis. A recent study from Italy revealed that the coexistence of ARGs and virulence factors could lead to the emergence of a new population of resistant clones, which poses a major challenge for addressing antimicrobial resistance ^[47]. We believe that P. mirabilis deserves more attention. In addition, CG3 has 77 shared ARG subtypes and 18 unique ARG subtypes. CG3 was widely distributed on all continents, and its major resistance genes were classified as β-lactams (e.g., cephalosporin and carbapenem), which is consistent with the prevalence of β-lactamase-carrying P. mirabilis strains reported in many countries and regions ^[48-50]. Therefore, attention should be paid to these CGs of P. mirabilis.

There are some limitations in this study that must be pointed out. First, our analysis of the ARGs and virulence genes of the CGs of P. mirabilis was only preliminary. We simply used the easy and fast bioinformatics pipeline (ABRicate) and did not perform relevant experimental studies such as gene validation using polymerase chain reaction and quantitative PCR. The main reason for this is that the sequences of all 695 P. mirabilis strains were downloaded from NCBI, which became the main reason why we were unable to perform many experiments. Second, the main aim of the cgMLST scheme is not only typing, but also outbreak tracing. We did not apply the cgMLST scheme to epidemiological practice in this study. Therefore, further research will be conducted in the future to confirm that our cgMLST scheme is well suited not only for molecular typing of P. mirabilis, but also for epidemiological practice of P. mirabilis. Despite these current limitations, this is still the first study that performed high-resolution molecular typing of P. mirabilis using WGS technology combined with bioinformatics pipelines (chewBBACA), our results support the prevention and control of P. mirabilis infections.

In summary, we carried out high-resolution molecular typing of P. mirabilis from all over the world and found that the CGs of P. mirabilis showed significant regional distribution differences. We identified a total of 205 CGs in 695 P. mirabilis strains. The proposed cgMLST scheme for P. mirabilis typed 100% (707/707) of the available isolates, and at least 95% of these selected loci were present in the genome. We also found that some P. mirabilis CGs carried a large number of shared ARGs (or virulence genotypes) and unique ARGs (or virulence genotypes). To the best of our knowledge, this is the first report of high-resolution molecular typing of P. mirabilis using whole genome sequencing technology combined with a bioinformatics pipeline (chewBBACA). We have published our scheme, associated codes, and files on GitHub (https://github.com/Natasha22222222/cgMLST-Proteus-mirabilis) to facilitate further discussion and contribute to the establishment of cgMLST for P. mirabilis.

The datasets for this study can be found in the NCBI (National Center for Biotechnology Information) genome sequence repository (http://www.ncbi.nlm.nih.gov/genome/). The RefSeq assembly accession of all included P. mirabilis sequences are provided in Supplementary Table S2.