Carbapenem-resistant P. rettgeri strain 18004577 was recovered from the urine sample of an 18-year-old male patient with chest and abdominal trauma who was admitted to Shandong Provincial Hospital, China, in 2018. The patient was diagnosed with multiple fractures and lung and kidney contusions after falling from a six-story building. Urinary tract infection occurred during hospitalization. Finally, the patient recovered and was discharged. The patient’s medication history and hospital course details were retrieved from the hospital record information system. The case history collection and reporting protocols used in the present study were approved by the Ethics Committee of Shandong Provincial Hospital.

Species identification was performed using matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (BioMérieux, France). Phenotypic detection of carbapenemases was conducted using the carbapenem inactivation method (CIM) and the EDTA-modified CIM (eCIM) test, according to the guidelines of the Clinical and Laboratory Standards Institute (CLSI). Antibiotic susceptibility testing for aztreonam, cefepime, ceftriaxone, ceftazidime, ertapenem, imipenem, piperacillin/tazobactam, trimethoprim/sulfamethoxazole, ciprofloxacin, levofloxacin, gentamicin, amikacin, ampicillin, ampicillin–sulbactam, cefazolin, cefotetan, tobramycin, and nitrofurantoin was conducted using the VITEK-2 compact system (BioMérieux, France); the test results were interpreted using CLSI breakpoints (CLSI, 2020).

Conjugation experiments were conducted as described previously, with sodium azide-resistant E.coli J53AziR as the recipient strain and P. rettgeri 18004577 as the donor strain. Transconjugants harboring carbapenemase resistance genes were selected on Mueller–Hinton agar plates containing 6 μg/ml ceftazidime and 100 mg/ml sodium azide. Antibiotic susceptibility tests and PCR were performed to confirm carbapenemase gene transfer (Poirel et al., 2011; Xiang et al., 2015).

S1-nuclease pulsed-field gel electrophoresis (S1-PFGE) was performed to determine the size and number of plasmids carried by P. rettgeri 18004577 (clinical strain) and J-18004577, a transconjugant. The genomic DNA of the clinical strain and transconjugant embedded in the gel plugs was digested with QuickCut S1 nuclease (Takara, Shiga, Japan) for 1 h or with QuickCut sfi nuclease (Takara, Shiga, Japan) for 2.5 h and then separated using S1-PFGE for 17 h. The pulse time was switched from 2.16 to 63.8 s. Salmonella strain H9812 was digested with QuickCut XbaI (Takara, Shiga, Japan) and used as the reference marker. Southern blotting was conducted to confirm the locations of blaNDM−1, blaOXA-10, and blaPER-4 using specific probes labeled with digoxigenin (Roche, Basel, Switzerland; Supplementary Figure 1A) A 621-bp probe for blaNDM-1 was synthesized using PCR amplification with the primers 5′-CGGAATGGCTCATCACGATCCAC-3′ (forward) and 5′-GGTTTGGCGATCTGGTTTTC-3′ (reverse). A 504-bp probe for blaOXA-10 was synthesized using the primers 5′-TCTGCCGAAGCCGTCAATGGT-3′ (forward) and 5′-ATATTCAGGTGCCGCCTCCGTTA-3′ (reverse). Finally, a 504-bp probe for blaPER-4 was synthesized using the primers 5`-GCAATACTCGGTCTCGCACAC′-3′ (forward) and 5′-TGATACGCAGTCTGAGCAACCT-3′ (reverse).

DNA was extracted from the clinical isolates and transformants using a genomic DNA commercial kit (Qiagen, Hilden, Germany). Genomic DNA was sequenced using the Illumina HiSeq platform (Novogene Co., Ltd., Beijing, China) and PacBio RSII sequencer (Biozeron Biological Technology Co., Ltd., Shanghai, China). Paired-end short Illumina reads were used to correct long PacBio reads using proofread for large-scale high-accuracy PacBio correction, and the corrected PacBio reads were then assembled de novo utilizing using the functions available at https://github.com/ruanjue/smartdenovo.

Sequence annotation was conducted using RAST 2.0¹ combined with BLASTP/BLASTN searches against UniProtKB/SwissProt and RefSeq databases. Annotation of the resistance genes and mobile elements was conducted using online databases, including CARD² and ISfinder.³ PHAge Search Tool Enhanced Release was used to identify prophages (Arndt et al., 2016). The virulence genes were identified by alignment of the gene sequences against the sequences available in the Virulence Factors Database (Liu et al., 2019).

The sequences of five plasmids were retrieved from the NCBI database for comparative analysis: pBML2531 (GenBank accession no. AP022376.1), pNDM15-1091 (GenBank accession no.CP012903.1), pRts1 (GenBank accession no. MN626604.1), pT-OXA-181 (GenBank accession no. NC_020123.1), and Rts1 (GenBank accession no. NC_003905.1).

A BLAST Ring Image Generator⁴ was used to generate and visualize the comparisons of the plasmids and their genetic structures. More detailed genome alignment between closely related plasmids was conducted using local BLAST, and the findings were visualized using Easyfig.⁵

Plasmid stability testing was performed using Luria–Bertani broth as previously described, but with some modifications (Nang et al., 2018). Briefly, the P. rettgeri 18004577 clinical strain was cultured at 37°C in a shaking bath (150 rpm/min) and serially passaged for 7 days with a 1:1,000 dilution of antibiotic-free Luria–Bertani broth. Then, 100 μl of the seventh-day culture was plated onto antibiotic-free Mueller Hinton agar. Subsequently, 50 colonies were randomly selected and subjected to PCR for amplification of the repA gene of the plasmid p18004577_Rts. The plasmids were considered stable if more than 85% of the colonies harbored the repA gene.

Combining with sequencing data of 29 published Providencia rettgeri genomes from the National Center for Biotechnology Information (NCBI) database, pan-genome analysis was carried out (Supplementary Table 1). The criteria for data-selection was that Providencia rettgeri strain with select columns “level” had “complete” genome files available on NCBI in the beginning of our project. We first determined the phylogenetic relationship between the 30 P. rettgeri strains using OrthoFinder (Emms and Kelly, 2015), by using the protein sequences of the strains obtained from the Prokaryotic Genome Annotation System (Prokka; Seemann, 2014). Further, the multiple sequence alignment analysis of the resulting single-copy direct-line homologous protein long sequence was carried out using the MAFFT software (Katoh, 2005). Then, we constructed the maximum-likelihood phylogenetic tree of 30 genomes using the RaxML-NG software, with a 1,000-bootstrap test (Kozlov et al., 2019). The calculation of the best amino-acid substitution model was performed using ModelTest-NG software before the phylogenetic tree construction (Darriba et al., 2020). After obtaining the best tree, we visualized the evolutionary tree using the R software with ggtree package (Yu, 2020). Average nucleotide identity (ANI) was performed to explore the taxonomic boundary of the genomes using FastANI (Jain et al., 2018).

After the determination of potential confounding strains, pan-genome analysis was conducted with Roary (Page et al., 2015) using the GFF3 files generated by Prokka (Seemann, 2014). To explore the molecular and biological functions of core genes from the 30 genomes, GO functional enrichment analysis was carried using the R language clusterProfiler package (Wu et al., 2021).

Carbapenem-resistant P. rettgeri 18004577 was identified as an MBL-producing strain using eCIM. The antimicrobial susceptibility testing results are presented in Table 1. The clinical strain was resistant to almost all of the 18 antibiotics, except for amikacin and trimethoprim/sulfamethoxazole. The transconjugant J-18004577 remained susceptible to amikacin, trimethoprim/sulfamethoxazole, aztreonam, and nitrofurantoin.

According to the whole-genome sequencing analysis, the complete genome of strain 18004577 contained a circular chromosome of 4.6 Mb with a G + C content of 41%, a circular plasmid p18004577_NDM of 273.2 Kb with a G + C content of 47.6%, and a circular plasmid p18004577_Rts of 146.2 Kb with a G + C content of 45.5%.

A total of 10 prophage regions were identified, of which eight were located on the genome and two were located on the plasmid p18004577_NDM (Supplementary Table 2). A CRISPR region was identified in the genome of 18004577 at nucleotide positions 565, 107–565, 213-bp, with one 29-bp long spacer sequence.

The distribution of virulence genes was investigated to identify the key pathogenicity genes of P. rettgeri 18004577. The strain harbored several virulence genes, such as flhA, fliC, clpB, hcp, fcl, and gmd, associated with flagellar biosynthesis, type VI secretion system, and O-antigen (Supplementary Table 3).

According to the results of S1-PFGE, the clinical strain P. rettgeri 18004577 harbored two plasmids of approximately 300 and 100 Kb (Supplementary Figure 1B). Based on the whole-genome analysis, the larger plasmid named p18004577_NDM (GenBank accession no. CP098041.1) was found to harbor multiple resistance genes, including bla_NDM-1, bla_OXA-10, bla_PER-4, aph(3′)-VI (kanamycin resistance), ant(2′′)-Ia (gentamicin and tobramycin resistance), ant(3′)-IIa (streptomycin resistance), sul1 (sulfamethoxazole resistance), catB8 and catA1 (chloramphenicol resistance), mph(E; erythromycin resistance), and tet (tetracycline resistance; Supplementary Table 4).

The New Delhi MBL (NDM)-harboring plasmid of the clinical strain was transferred into E. coli J53Azi^R via conjugation. The presence of NDM-1 in the transconjugant was confirmed using PCR. Two resistance genes were identified in the transconjugant J-18004577––bla_NDM-1 and aph(3′)-VI––which were obtained from the clinical isolate P. rettgeri 18004577 by conjugation assay (Figure 1). Southern blotting confirmed that bla_NDM-1, bla_OXA-10, and bla_PER-4 were located on the larger plasmid in the clinical strain, and only bla_NDM-1 was present on the hybrid plasmid, with a size of approximately 220 Kb of the transconjugant (Supplementary Figures 1C,D).

p18004577_Rts was found to be a type IncT plasmid with a 146,248-bp size and an average GC content of 45.5%; it had at least 202 predicted coding sequences. No resistance gene was predicted, and it included 15 tra-type genes, except for trb, in the transfer region. BLAST analysis showed that p18004577_Rts shared 99% nucleotide identity with 98% cover with plasmid Rts1 from Proteus vulgaris (NC_003905.1; Figure 2).

p18004577_Rts (GenBank accession no. CP098042.1) harbored by P. rettgeri 18004577 contained a replication initiation protein encoded by repA (at 155–1177 bp) and a short segment containing the replication origin ori (at 51–238 bp) that provides two elements for its autonomous replication (Murata et al., 2002), a process which is similar to the one found in plasmid Rts1. The replication of this mini-Rts1 plasmid was stable at 37°C (Itoh et al., 1982) The stability of p18004577_Rts was confirmed using PCR amplification performed after 7 days of culture; 48 of 50 colonies were found to carry repA, indicating the high stability of p18004577_Rts.

pJ18004577_NDM (GenBank accession no. CP114206) was a hybrid plasmid generated from p18004577_NDM to p18004577_Rts, with a length of 154,303 bp, a GC content of 45.5%, and 209 predicted coding sequences. The features of the variable region (at 8,779–24,607 bp) of pJ18004577_NDM shared high homology with the bla_NDM-1-containing region of p18004577_NDM, and the backbone region of pJ18004577_NDM was almost identical to that of p18004577_Rts (Figure 3).

According to the whole-genome analysis, p18004577_NDM harboring bla_NDM-1 from clinical strain P. rettgeri 18004577 was a 273,271-bp long untypeable plasmid, based on replication module analysis. p18004577_NDM contained an accessory multidrug-resistant region generated by ΔTn1696, class 1 integrons (IntI) harboring bla_OXA-10, ant (2′)-Ia, catB8, bla_PER-4 and sul, ΔIS1326, ΔTn125 harboring bla_NDM-1, ΔISKoX2, ΔTn21, mer operon, and ΔTn2501, which conferred resistance to aminoglycosides, quinolones, and β-lactams. In addition, two prophage regions (at 218,975–233,789 and 288,401–260,814 bp) were predicted, the physiological functions of which are unknown (Figure 4).

In p18004577_NDM, a variant of Tn1696 was observed upstream of the class I integron, which was found to be similar to pSx1 (accession no. CP013115), pT_OXA_181 (accession no. JQ996150), and Tn1696 (accession no. U12338; Figure 5). The class I integron (at 27,283–48,006 bp) was flanked by an IS26 sequence that interrupted the transposase TniA at a downstream position, and the cassette content of the class I integron comprised multiple resistance genes, such as bla_OXA-10, ant (2′)-Ia, catB8, bla_PER-4, two Δsul1, two qacEdelta1, and one IS1326 element truncated by IS1353. The bla_NDM-1 containing region (at 53,799–58,491 bp) was truncated by Tn125 that was flanked by an IS14 element at an upstream location and by ISKoX2 inserted by IS26 at a downstream location (Figures 1, 4).

The plasmid pJ18004577_NDM showed partial homology with several other plasmids, such as pBML2531, pNDM15-1,091, pRts1, pT-OXA-181, and Rts1 (Figure 3). The features of the variable region (at 8,779–24,607 bp) of pJ18004577_NDM shared high homology with the bla_NDM-1-containing region of p18004577_NDM, and the backbone region of pJ18004577_NDM was highly similar to that of p18004577_Rts.

Based on the results of genetic comparison, the main possible formation process of pJ18004577_NDM was the insertion of the [ΔISKox2-IS26-ΔISKox2]-[aph(3′)-VI]-bla_NDM-1 module from p18004577_NDM (at 48,007–63,004 bp) into plasmid p18004577_Rts (at 8,779–17,992 bp; Figure 1). The insertion module in pJ18004577_NDM seemed to be surrounded by an 8-bp target site duplication with the sequence “GCTGAGAT” (at 8,871–8,878 bp) and “TTGAAGAA” (at 13,292–13,299 bp) at random sites; this may have been due to a cointegrate (Figure 6).

OrthoFinder assigned 123,540 genes (98.3% of the total) to 6,649 orthogroups. Fifty percent of all genes were in orthogroups with 30 or more genes (G50 was 30) and were contained in the largest 2042 orthogroups (O50 was 2042). There were 2,591 orthogroups with all species present and 2,442 of these consisted entirely of single-copy genes.

The maximum-likelihood phylogenetic tree constructed based on the best amino acid substitution model of 30 genomes single-copy lineal homolog proteins was used to further explore genomic similarities among strains showed in Figure 7A. The analysis of 30 P. rettgeri genomes by FastANI showed that higher the association between the microbial species, the higher the ANI values between each other (Figure 7B). The ML phylogenetic tree revealed that strains could be broadly clustered into two major clades base on the genetic distance calculated by the ANI value. The sample sources were scattered among the two clades of phylogenetic trees.

Roary analysis of P. rettgeri pan-genome identified 415 core, 756 soft core, 5,744 shell and 12,967 cloud genes. As the number of sequenced genomes increases, the species’ pan-genome size converges to a certain value, highlighting the “close” nature of P. rettgeri pan-genome. Accordingly, we obtained four different classes of genes belonging to “core” (29 ≤ strains ≤ 30), “soft core” (28 ≤ strains < 29), “shell” (4 ≤ strains < 28), and “cloud” (strains <4) groups, respectively (Figure 8A). Pan-genomics of the whole dataset was performed, followed by plotting the pan genome matrix. The matrix disclosed the deviance in the presence-absence profile of the 30 P. rettgeri strains (Figures 8B–D). The functional enrichment analysis of these 415 core genomes with the results are shown in Supplementary Figure 2.