Genomic Context of the Two Integrons of ST-111 Pseudomonas aeruginosa AG1: A VIM-2-carrying Old-acquaintance and a Novel Imp-18-carrying Integron

Pseudomonas aeruginosa is an opportunist and versatile organism responsible for infections mainly in immunocompromised hosts. This pathogen has high intrinsic resistance to most antimicrobials. P. aeruginosa AG1 (PaeAG1) is a Costa Rican high-risk ST-111 strain with resistance to multiple antibiotics, including carbapenems, due to the activity of VIM-2 and IMP-18 metallo-β-lactamases (MBLs). These genes are harbored in two class 1 integrons located inone out of the 57 PaeAG1 genomic islands. However, the genomic context associated to these determinants in PaeAG1 and other P. aeruginosa strains is unclear. Thus, we rst assessed the transcriptional activity of VIM-2 and IMP-18 genes when exposed to imipenem (a carbapenem) by RT-qPCR. To select related genomes to PaeAG1, we implemented a pan-genome analysis to dene and up-date the phylogenetic relationship among complete P. aeruginosa genomes. We also studied the PaeAG1 genomic islands content in the related strains and nally we described the architecture and possible evolutionary steps of the genomic regions around the VIM-2- and IMP-18-carrying integrons. associated with the resistance to carbapenems in this high-risk P. aeruginosa using comparative genomics.


Introduction
Pseudomonas aeruginosa is an opportunist and versatile pathogen able to survive in a wide variety of environments (Klockgether et al., 2010). With a large genome (6-7.5 Mb), a considerable proportion of the genome (>8%) of P. aeruginosa strains is dedicated to regulatory functions (Cabot et al., 2016) resulting in a consequent diversity of metabolic capabilities and responses to stress (J. A. . Due to these characteristics P. aeruginosa is responsible for infections among immunocompromised hosts (Lu et al., 2016) and nosocomial infections (Fernández, Corral-Lugo, & Krell, 2018). This pathogen has high intrinsic resistance to most clinically used antimicrobials (Brazas, Brazas, Hancock, & Hancock, 2005), many of them by multidrug-resistant or extensively drug-resistant strains (Oliver, Mulet, López-Causapé, & Juan, 2015). This severely compromises the selection of appropriate treatments (X.  leading to a signi cant increase in the morbidity and mortality of these infections. According to World Health Organization (WHO), resistance to carbapenems in P. aeruginosa, Acinetobacter baumannii, and Enterobacteriaceae family is considered a critical issue in the context of antibiotic resistance, being classi ed as Priority 1 group (World Health Organization, 2017). P. aeruginosa AG1 (PaeAG1) is a particular P. aeruginosa strain isolated from an immunocompromised patient in a Costa Rican hospital in 2010 (Toval et al., 2015). This strain is resistant to multiple antibiotics such as β-lactams (including carbapenems), aminoglycosides, and uoroquinolones, being only sensitive to colistin. This was the rst reported strain of a P. aeruginosa isolate carrying both VIM-2 and IMP-18 genes encoding for metallo-β-lactamases (MBLs) enzymes, both with carbapenemase activity (Toval et al., 2015). As demonstrated in our previous works, including the genome assembly (GenBank CP045739) (J.-A. Molina-Mora, Campos-Sánchez, Rodríguez, Shi, & García, 2020), these genes belong to two independent class 1 integrons, each contained in one out of the 57 predicted genomic islands of PaeAG1 (J.-A. Toval et al., 2015). Other elements such as six phages, mobile genetic elements, and some virulence factors are also harbored in genomic islands. Cipro oxacin exposure in PaeAG1 induces phage activity with a very complex activity, affecting the growth in spite of this strain being resistant to this antibiotic (J. A. . In addition, PaeAG1 has a not functional CRISPR-Cas system and molecular genotyping by multilocus sequence type (MLST) classi es PaeAG1 as a high-risk sequence type 111 (ST-111) strain (J.-A. . ST-111 is a lineage that belongs to the O12 serotype, including a multi-resistance pro le and the ability to colonize nosocomial environments (X. Turton et al., 2015;Witney et al., 2014;Woodford, Turton, & Livermore, 2011). Toghether with ST-235 and ST-175 genotypes, ST-111 belongs to the high-risk group in P. aeruginosa (Oliver et al., 2015). High-risk clones are frequently associated with epidemics where multidrug resistance impairs thetreatment (Petitjean et al., 2017).
In this context, it is considered that P. aeruginosa high-risk clones are part of a non-clonal epidemic population structure (Oliver et al., 2015;Petitjean et al., 2017), many carrying genomic determinants such as carbapenemases or extended-spectrum β-lactamases (Oliver et al., 2015). Carbapenemases include Ambler class A enzymes such as KPC and GES variants, Ambler class B MBLs (IMP, VIM, SPM, GIM, NDM, and FIM type), and Amber class D (OXA variants) enzymes (Farajzadeh Sheikh et al., 2019;Hong et al., 2015). In Costa Rica, isolation of carbapenem-resistant P. aeruginosa strains is relatively common in some major hospitals as we reported, most of them carrying VIM or IMP alleles and up to a 63.1% prevalence (Toval et al., 2015). This is much higher than the frequencies observed in other countries (Hong et al., 2015).
VIM and IMP genes, as well as other MBLs, are frequently found as part of gene cassettes carried in integrons (Walsh, 2005;Zhao & Hu, 2011), leading to the dissemination of multidrug resistance among Gram-negative bacteria (Jones-Dias et al., 2016). Thus, there is a growing interest in the reconstruction of the genomic context of mobile elements (in particular for integrons) to gain insights into bacterial evolution and its association with human activities, as well as to identify possible ways to mitigate antibiotic resistance (Ghaly, Chow, Asher, Waldron, & Gillings, 2017). However, there are few studies about the genomic context of P. aeruginosa high-risk clones associated with integrons (Chowdhury et al., 2016).
In this sense, comparative genomic strategies can provide insights not only about gene content, architecture and evolutionary details, but also dynamics of mobile genetic elements, pathogenicity determinants, and others (Peter et al., 2019). Several studies at the genomic level have been implemented to describe the molecular diversity in P. aeruginosa (including high-risk clones) using different comparative approaches (Xavier Petitjean et al., 2017;Turton et al., 2015). In our previous work, a pan-genome comparison was performed using only 11 ST-111 strains, in which the gene content was used to phylogenetically describe PaeAG1. However, as discussed there, an in-depth analysis including a signi cant number of strains (not only 11) is required to better understand the genomic features of PaeAG1 and other ST-111 strains (J.-A. Molina-Mora et al., 2020).
Since PaeAG1 has special genomic features regarding antibiotic multi-resistance, including VIM-2 and IMP-18 genes with carbapenemase activity, 57 genomic islands, and a ST-111 pro le, we hypothesized that the comparative genomics can reveal insights about the evolution and landscape of genomic regions around the MBLs-carrying integrons of PaeAG1. Thus, the aim of the study was to compare the PaeAG1 genome against other P. aeruginosa sequences using comparative genomics to describe phylogenetic relationships, genomic islands content, and architecture of genomic regions associated with the VIM-2-and IMP-18-carrying integrons of PaeAG1. We rst demonstrated that VIM-2 and IMP-18 are functional genes that can be induced after treatment with imipenem (a carbapenem antibiotic). We then analyzed all the complete P. aeruginosa genomes using a pan-genome analysis approach to identify related genomes to PaeAG1, revealing that whole genome sequences can separate clones by MLST pro le (ST). Afterward, PaeAG1 genomic islands were searched in the related genomes, including all the ST-111 genomes, and diverse presence/absence patterns were found in related genomes. Finally, speci c genomic regions associated with the two integrons were reconstructed and characterized to compare the gene content and architecture in close genomes. The genomic region associated with the VIM-2-carrying integron (In59-like) was completely found in other two ST-111 strains (i.e. it is an old-acquaintance integron), but an IMP-18-carrying integron (registered as In1666), with an architecture never reported before, was found when the landscape of the related genomic region was described.

Bacterial isolate
The PaeAG1 strain is a Costa Rican isolate with resistance to β-lactams (including carbapenems, MIC Imipenem >32 µg/mL), aminoglycosides, and uoroquinolones, being only sensitive to colistin. We recently assembled and annotated the PaeAG1 genome (J.-A.  and data is available in Genbank under accession CP045739 (Bioproject PRJNA587210).
2.2 RT-qPCR for VIM-2 and IMP-18 expression after imipenem exposure In order to study the expression of VIM-2 and IMP-18 genes by imipenem exposure in PaeAG1, experiments of growth curves and RT-qPCR were performed.
Growth curves assay: Three aliquots of pre-cultured PaeAG1 cells were added to fresh Lysogenic Broth (LB) broth to an initial optical density (OD 600nm ) of 0.01. Each aliquot was treated with 0.0 (control), 25.0 or 50.0 µg/mL of imipenem. Growth was monitored at times 0, 2, 4,6,8,12, and 16 hours. The assay was performed in triplicates. Two speci c aliquots at times 6 and 12 hours were taken for RT-qPCR assay, as follows.
RNA isolation: Aliquots at times 6 and 12 hours after imipenem exposure were preserved using the RNA protect reagent (QIAGEN). Total RNA was extracted using the RNeasy Mini kit (QIAGEN, UK) following the manufacturer´s instructions. Subsequently, RNA was transcribed into cDNA with the Maxima H Minus First Strand cDNA Synthesis kit (Thermo Scienti c™ Inc.). In the different steps, the quality and quantity of extracted RNA or cDNA were determined using a Nanodrop (Nanodrop 2000, Thermo Scienti c™ Inc.).
RT-qPCR: The standard curve method was implemented to quantify the expression of target and reference genes. Each reaction mixture contained 12.5 μL of SYBR green Master Mix (Thermo Scienti c™ Inc.), 0.25 μL of each primer, 10 μL of PCR-grade water, and 2 μL of cDNA. Thermocycling was performed in a StepOnePlus Real-Time PCR System (Thermo Scienti c™ Inc.). For VIM-2 and IMP-18 genes, the assay was run with denaturation at 95°C (10 min), 35 ampli cation cycles of 94°C (20 s), 53°C (45 s), and 72°C (30 s), with data acquisition at 72°C. For the rpoD gene, conditions were denaturation at 95°C (10 min), 45 ampli cation cycles of 95°C (15 s), 20°C (10 s), and 72°C (15 s), with data acquisition at 72°C. Melt curve data were used to determine the identity of the amplicons. .
Relative gene expression analysis: Gene expression of VIM-2 and IMP-18 in the experimental conditions (0, 25, and 50 μg/mL imipenem) were normalized using the rpoD housekeeping gene. Data were analyzed using the delta-delta Ct method (12). The change in gene expression within samples (time and antibiotic concentration) was calculated with reference to the control (0 μg/mL imipenem) and a two-way ANOVA test was performed between conditions (95% con dence level).

Datasets of complete P. aeruginosa genome sequences
To compare all the complete genomic sequences of P. aeruginosa by a pan-genome analysis, metadata (including strain names, alternative ID, gene content, MLST pro le, and others), genome and protein sequences (".fasta" format), and annotation (Genbank ".gbk" and ".tab" formats) les were retrieved from Pseudomonas Genomes Database (PGDB, https://pseudomonas.com).

Comparative genomic analysis by a pan-genome approach
Since differences in annotation were identi ed for many sequences, even in the same genomic regions, we decided to identify and annotate genes from the complete genomic sequences using the same approach. To achieve this, gene prediction and annotation were done using Prokka v1.13.3 (with --genus Pseudomonas --species aeruginosa and other parameters by default con guration) (Seemann, 2014). The Prokka annotation les (in ".gbk" format) were used to run the phylogenetic analysis by a pan-genome approach based on gene content in the Roary program v3.12.0 (Page et al., 2015) with default parameters. The phylogenetic tree (".newick" le) was visualized using Interactive Tree Of Life Tool (iTOL, https://itol.embl.de/) v5 (Letunic & Bork, 2019), and strain names and MLST pro les were incorporated for each strain. For strains with unknown MLST, the pro le was veri ed using the complete genome sequence approach (Larsen et al., 2012) in the MLST tool v2.0 (https://cge.cbs.dtu.dk/services/MLST/). For a functional analysis for all core-genes, STRINGdb (https://string-db.org/) was used to identify signi cantly enriched KEGG pathways (cutoff of false discovery rate FDR < 0.05).

Comparative analysis of the presence of PaeAG1 genomic islands in other strains
The 57 PaeAG1 genomic islands were previously identi ed using IslandViewer v4 (www.pathogenomics.sfu.ca/islandviewer/), as we reported recently (J.-A. . Regions of the genomic islands (".bed" le) were downloaded from the same platform and sequences (".fasta" format) were obtained using the getfasta function in bedtools software v2.29.2 (Quinlan & Hall, 2010). The distribution of genomic islands along the genome was visualized using the BLAST Ring Image Generator BRIG tool v0.95 (Alikhan, Petty, Ben Zakour, & Beatson, 2011).
In order to determine the presence and frequency of these genomic islands in other strains, a comparative analysis based on sequence alignment was done. Thus, we implemented a BLASTn pipeline to align PaeAG1 genomic island sequences and the complete genome sequences of all strains. A minimum length for coverage of 95% (overlap between query and subject sequences) and 80% of minimum sequence identity between sequences were used to de ne that a speci c genomic island was present in a strain, similar to other studies (Holmes, Dallman, Shabaan, Hanson, & Allison, 2018;Kluytmans-van den Bergh et al., 2016). A nal comparison of the presence/absence of genomic islands was done for selected strains (see Results) using a small phylogenetic tree and a heatmap, which were visualized using phylo.heatmap function from phytools package v0.7-20 (https://www.rdocumentation.org/packages/phytools), in the R software (https://www.r-project.org/).

Landscape of genomic regions associated with the two class 1 integrons of PaeAG1
Two complete and independent class 1 integrons were previously identi ed in PaeAG1, one carrying the VIM-2 gene and another harboring the IMP-18 gene (J.-A. Molina-Mora et al., 2020). To better understand the possible evolutionary history of these integrons and their potential for lateral transfer, we reconstructed the genetic landscape of the genomic regions around these elements. The identity of the integrons was investigated using INTEGRALL database (http://integrall.bio.ua.pt). For the new integron (see Results), the same database was used for the registry and the integron number assignment.
Since the two integrons are absent in the reference strain Pae-PAO1, an alignment of the genomic regions (BLASTn) and another of amino acid (AA) sequences (BLASTp) were used to identify the limits of the complete inserted region in PaeAG1. The two speci c inserted regions were composed of two or more genomic islands in a row, as obtained in our previous study (grouped or with overlapping regions) (J.-A. Molina-Mora et al., 2020). Thus, regions were called GIC VIM-2 (genomic island cluster containing VIM-2carrying integron) and GIC IMP-18 (genomic island cluster harboring the IMP-18-carrying integron).
A nal alignment (BLASTn) of the extended regions (three coding genes on each side) of GIC VIM-2 and GIC IMP-18 was done against selected genomes. Genomes selection was done based on the phylogenetic relationships of strains close to PaeAG1 (pan-genome analysis) and the pro le of presence/absence of the PaeAG1 genomic islands in other strains. All the syntenic regions of selected strains were compared using annotation les in Easy g software v2.2.3 (Sullivan, Petty, & Beatson, 2011), leading to visualize alignments, gene content and identity, exclusive/shared elements by strain and possible evolutionary steps, and others.

Expression of VIM-2 and IMP-18 genes is induced after imipenem treatment in PaeAG1
In order to assess the functional activity of VIM-2 and IMP-18 genes, a RT-qPCR was performed. Exposure to imipenem had no effects on the growth curves of PaeAG1 ( Fig. 1-A). Evaluation of gene expression after exposure to imipenem ( Fig. 1-B-C) showed that VIM-2 and IMP-18 increased its expression at least by a 1.7-fold (with respect to control) at 6 hours, but only 1.1-fold at 12 hours. This observation was independent of the imipenem concentration (25 or 50 μg/mL), as supported by the statistical analysis in which changes in the relative expression by time but not by concentration were signi cant for each gene.
3.2 Pan-genome analysis with the complete genome sequences de nes P. aeruginosa clusters which correlates with the MLST genotyping pro le To select related genomes to PaeAG1, a total of 211 strains were analyzed to compare the genomic composition including PaeAG1 (Supplementary le 1 All_strains_information.xlsx contains the list of all the selected genomes, ID, strain, MLST pro le, and others). Gene content comparison was done based on a pan-genome approach. A total of 2726 genes were identi ed as part of the core-genome (present > 99% strains). More details of results and complementary plots are provided in the Supplementary le 2 Pangenome analysis results.xlsx. Enrichment analysis of KEGG pathways for all core genes (Table S1) found 42 biological processes implicated in several metabolism routes related to energy (carbon, fatty acids, amino acids), DNA and RNA, ribosomal activity, protein synthesis, and others.
As shown in Fig. 2, the similarity in the genomic composition by pan-genome analysis de nes a phylogenetic tree able to separate groups that can be described in turn by the MLST genotyping pro le. Although we identi ed a total of 67 different MLST pro les (and unknown cases), many of them resulted in low frequency. For example, 35 different ST classes had only a single strain (35 strains, 17% of all genomes) and 88 strains (42%) belonged to the 56 ST pro les with less than ve genomes. Also, 44 strains (21%) had an allelic composition with an unknown ST pro le. On the other hand, a total of 79 (37%) genomes corresponded to 11 ST classes with ve or more strains. The last was evidenced using different colors by ST pro le (as shown in Fig. 2), meanwhile strains belonging to low-frequency ST pro les were colored in the same way. Representative genomes such as the reference strain Pae-PAO1 (ST-549, purple cluster) and Pae-UCBPP-PA14 (ST-253, yellow group) were identi ed in the main ST groups.
PaeAG1 was located in the same group with the other nine ST-111 strains in a separated cluster (green). Other two ST pro les (low-frequency ST-234 and ST-654) and one unknown case (Pae-Pa84 strain) were kept close to this group. The whole group of these related strains, and the reference strain Pae-PAO1, were used for subsequent analysis, including their phylogenetic relationships. For other high-risk clones, a single ST-175 genome was identi ed, and a clear cluster was found for the ten ST-235 genomes (including other genomes with unknown pro les).
3.3 Varying pro les of the presence/absence of the 57 PaeAG1 genomic islands are found in the ST-111 strains and related genomes A comparative analysis based on sequence alignment was run to determine the presence and frequency of the PaeAG1 genomic islands in other phylogenetically related strains. Genomic islands locus were previously predicted (J.-A. Molina-Mora et al., 2020). We rst represented the distribution of the genomic islands along the PaeAG1 genome, as presented in Fig. 3. Many of the islands are kept together, including overlapping regions or an arrangement in a row. Thus, we termed this as a genomic islands cluster (GIC) to refer to this group of islands. In Fig. 3, GICs correspond to the genomic regions labeled as joined names of the genomic islands, for example "GI48-49" represents the genomic region of islands GI48 and GI49. In some cases, each genomic island in the cluster can be differentially distributed in the genomes (for example GI48 is present in PaeAG1 and Pae-97, but GI49 is only found in PaeAG1, Fig. 4). For this reason, we do not re-de ne the locus neither joined the islands.
Analysis of the presence/absence of PaeAG1 genomic islands in other ST-111 strains and related genomes is shown in Fig. 4. Pro les for all the 211 is available in the Supplementary le 1 All_strains_information.xlsx, including total counts of strains by genomic islands, and total genomic islands per genome. The closest genomes to PaeAG1 (Pae-RIVM-EMC2982 and Pae-Carb0163) had the most similar pro les in the genomic islands content (41 genomic islands), but different patterns are obtained for other ST-111 strains. None of the islands is present in the reference genome Pae-PAO1, and other few genomic islands are rarely present in other non ST-111 strains.
On the other hand, two particular genomic islands were particularly recognized due to they carry the two PaeAG1 integrons. GI27 genomic island harbors the VIM-2-carrying integron, while IMP-18-carrying integron belongs to GI49. As shown in Fig. 4, GI27 (red) is present in PaeAG1 and two other ST-111 strains, and it is also absent in the rest of the 208 genomes. GI49 (blue) is unique to PaeAG1 and it is not it is present in none of the other 210 strains in the study.
Based on phylogenetic relationships, ST pro le, and genomic islands content, we selected speci c genomes to compare the GICs associated with the integrons. As shown in Fig. 4, the four genomic islands of GIC VIM-2 (GI27-30) are differentially present in the genomes. For example, GI28 and G29 are present in eight strains, but GI27 in three and G30 in four. To speci cally compare the genomic regions of GIC VIM-2 , we used the reference Pae-PAO1, Pae-RIVM-EMC2982 (with the four genomic islands), and Pae-AR445 (with three of the genomic islands). For the case of GIC IMP-18 , the two islands GI48 and GI49 are absent in other ST-111 strains, but GI48 is present in Pae-97. Except for this case, no other strains in all 211 genomes were identi ed harboring both islands. To compare the genomic regions, the reference genome Pae-PAO1, Pae-RIVM-EMC2982 as the closest genome, and Pae-97 (the only genome sharing a section of the GIC) were used.

GIC VIM-2 is a known region containing the old-acquaintance VIM-2-carrying integron in PaeAG1
With the aim of describing the possible evolutionary history of the VIM-2-carrying integron in PaeAG1, we described the architecture of the genomic regions delimited by the GIC VIM-2 (including three genes on each side: 35 798 bp and 32 protein-coding genes). Using Pae-PAO1 as reference, we found that genomic insertion occurred in the middle of the PA2229 gene, as shown in the top of Fig. 5. The insertion resulted mostly present in Pae-AR445 (coverage 94% and identity 99.97% of the PaeAG1 region), but without most of the integron (integrase intI1 and sul1 are present, unlike the gene cassette including VIM-2). However, a full coverage region was identi ed in Pae-RIV-EMC2982, with a 100% coverage and identity 99.99%. The only two variants identi ed in the full region were non-synonymous mutations, with an amino-acid change in PaeAG1_03254 (transcriptional regulator merD, 99.0% identity) and PaeAG1_03255 (mercuric reductase merA, 99.8% identity). See Table 2 and Supplementary Table S2 for more details. Although not shown in Fig. 5, alignment was also done for Pae-Carb0163, which has the same pro le of genomic islands content as Pae-RIV-EMC2982. In this case, a 100% coverage and identity 99.87% (45 variants) were obtained in the GIC VIM-2 region; most of the variants resulted in a change in the amino-acid sequence in PaeAG1_03245 (aacA29a, part of the integron with a 95.8% identity resulting in aacA29e allele), but also affecting other three proteins (mercuric reductase, integrase IntI, and a transposase). See Supplementary Table S2 for more details.
Regarding the gene content (Table 2), this genomic insertion contains the complete integron carrying VIM-2 gene. The composition of this integron is described in Fig. 5 (bottom), containing classical elements int1, attI, sul1, and the gene cassette (with aacA29a-b and VIM-2) of a class 1 integron, being classi ed as In59-like. Furthermore, GIC VIM-2 has at least other mobile genetic elements, including transposases and recombinases modules. Other coding modules are associated with mercury metabolism or they remain unknown (hypothetical proteins). Details of the protein alignment of PaeAG1 against four genomes are also provided (Supplementary Table S2). Reconstruction of the evolutionary steps related to the conformation of this genomic region includes the participation of four transposons (Tn402, Tn21-like, a disrupted and another complete Tn4661) as shown in Fig. 7-A. See details in Discussion.
Considering the full coverage and very high identity in at least two genomes, Pae-RIVM-EMC2982 and Pae-Carb0163, GIC VIM-2 can be considered a genomic region present in two well-known VIM-2+ strains, being this gene located in an old-acquaintance class 1 integron (In59-like).
In the case of other VIM-2-carrying integrons, more than 90 elements have been reported (according to INTEGRALL database) and most of them include aacA alleles, as in the case of the In59-like integron present in PaeAG1 and other strains. A similar arrangement of this element is reported for other integrons such as In103, In927, In1025, and In1163. Details of all the 93 VIM-2-carrying integrons are described in Supplementary Table S3. 3.5 GIC IMP-18 is a PaeAG1 exclusive genomic region harboring a new IMP-18-carrying integron Similarly, we compared four genomes to describe the architecture of the genomic regions delimited by the GIC IMP-18 (including three genes on each side: 30 258 bp and 29 protein-coding genes). Using Pae-PAO1 as reference, we found that genomic insertion occurred between the genes PA4704 and PA4705, as shown in the top of Fig. 6. Genomic islands GI48-49 are absent in Pae-RIV-EMC2982 and Pae-Carb0163 genomes.
BLAST of GIC IMP-18 identi ed the highest scored sequence in the Pae-97 genome (ST-234). Thus, since Pae-97 carries GI48, the syntenic comparison was done using this genome (Fig. 6). The analysis revealed a 77% coverage with identity 99.92%. The Pae-97 integron also contains Int1, aacA genes, and another allele of the IMP gene (IMP-1), all with a different arrangement.
Regarding gene content (Table 3), this genomic insertion contains the complete integron carrying the IMP-18 gene. The composition of this integron is described in Fig. 6 (bottom), containing int1, attI, sul1, and the gene cassette (IMP-18, gcuD, OXA-2, and aacA4). GIC IMP-18 also has genes coding for endonucleases and recombinases, or hypothetical proteins. Details of the protein alignment of PaeAG1 against the four genomes are also provided (see Supplementary Table S4).
Considering the absence of the complete region in other genomes and the rst report of the architecture of this integron, GIC IMP-18 can be considered a PaeAG1 exclusive region harboring a new IMP-18-carrying integron. This integron was registered as In1666 in INTEGRALL database. According to the same database, at least seven IMP-18-carrying integrons have been reported. OXA, gcuD, aacA and/or aadA genes are frequently present with different arrangements, as shown in Table 4, including the new integron In1666.
Conformation of GIC IMP-18 region seems to include the participation of at least three mobile elements (the new integron In1666, insertion sequence IS1326, and transposon TnAs3) as shown in Fig. 7-B. However, a lack of information about the role of other elements (regions without matching sequences) makes it di cult to complete the possible evolutionary steps related to this genomic region.
In summary, the pan-genome analysis leads us to identify that the genomic content can separate groups according to the ST pro le (MLST genotyping). All the ST-111 strains, including PaeAG1, resulted in the same phylogenetic group but different presence/absence pro les of PaeAG1 genomic islands were identi ed in other strains, even for grouped genomic islands, the GICs. Analysis of the landscape of regions GIC VIM-2 and GIC IMP-18 revealed one known and another new arrangement of genomic sequences in PaeAG1, harboring two independent MBLs-carrying integrons. The IMP-18-carrying integron has a unique and exclusive composition, reported here for the rst time.

Discussion
Multi-resistance to antibiotics is a major threat to the public health because of the continuous emergence, worldwide spread, and increas ing prevalence (Hong et al., 2015). With a high-risk ST-111 pro le, PaeAG1 is a critical organism due to its resistance to multiple antibiotics, in particular to carbapenems (World Health Organization, 2017). In our previous work, the pan-genome analysis was focused on the gene content using only 11 genomes (ST-111) to phylogenetically compare PaeAG1 and other ST-111 strains (J.-A. Molina-Mora et al., 2020). In the current work, we included all the P. aeruginosa strains with complete genomes (211 strains) to update the core genome for this group, as well as we characterized the genomic islands and the genomic context of the two integrons, including the possible evolutionary steps that can explain the current state of the genomic architecture. Thus, these new results give more insights into the biology and evolution of the ST-111 strains, in particular for PaeAG1.
We rst demonstrated that the expression of VIM-2 and IMP-18 genes (with carbapenemase activity) is induced after imipenem exposure, evidencing that are functional genes. A similar pattern in the expression levels, with an increment at 6 hours and a reduction at 12 hours, was found for both genes after treatment. Scarce studies have studied the expression patterns of various MBLs at the same time, including one recent study which found similar expression between VIM, IMP, and other MBL enzymes but depending on the sensitivity pro le (Singh et al., 2019). Although the similar expression levels of VIM-2 and IMP-18 could be obtained by chance, a mechanistic relationship cannot be discarded. For example, an imipenem effect on PaeAG1 could be related to the induction of integrons as a response to stress. Cell stress after exposure to antibiotics is known to induce the SOS response, and the gene cassettes of the integrons can subsequently be induced (Hocquet et al., 2012). Since the two integrons of PaeAG1 are class 1 elements, they could be regulated similarly. More analyses are needed to validate or discard this possibility.
In the pan-genome analysis, we were able not only to reveal that 211 whole-genome sequences could separate clones by ST pro le (MLST) but also the identi cation of core and accessory genes was achieved. Other pan-genome analyses in P. aeruginosa also found clusters that could be identi ed by the ST pro le (Aguilar-Rodea et al., 2017;Weiser et al., 2019). While multiple comparative genomic analyses (many using a pan-genome approach) have been reported for P. aeruginosa (Aguilar-Rodea et al., 2017;Chowdhury et al., 2016;Freschi et al., 2019;Gomila, Peña, Mulet, Lalucat, & García-Valdés, 2015;Hilker et al., 2015;Mosquera-Rendón et al., 2016;Ozer, Allen, & Hauser, 2014;Poulsen et al., 2019;Valot et al., 2015;Weiser et al., 2019;Wendt & Heo, 2016), most of them include incomplete, fragmented or draft genomes, or sequences of few genes. In 2015, complete genomes were used in a similar approach, but only 17 genomes were available (NCBI), and only three corresponded to high-risk clones (Valot et al., 2015). Thus, our analysis provides an update of the general status of relationships of the 211 available complete genomes by pan-genome analysis.
In our approach, only complete genome sequences were used, contrasting with most of the other pangenome studies. Thus, the analyses were based on high con dent sequences, which can be considered a strong point of our strategy. The previous study using complete genomes included the only 17 available sequences in 2015 (Valot et al., 2015). In the same line, another advantage of our approach is that we used the whole gene content, unlike other studies that have used few genes (Aguilar-Rodea et al., 2017) or variant analyses (Weiser et al., 2019) to relate strains. The last is a particular consideration for ubiquitous bacteria with a very plastic genome. A limitation of our approach is that the number of complete genomes is relatively small in comparison to the available draft genomes (>4500 according to www.pseudomonas.com). Well-assembled draft genomes could be considered and more information could be generated, but a robust selection must be applied. However, that option was not considered here.
Biologically, the relatively high number of conserved genes in the core-genome can be associated with the ability to conquer multiple environments and to facilitate infectious capability towards a large set of hosts (Valot et al., 2015). According to functional analysis, 42 KEGG pathways (energy metabolism, nucleic acids, amino acids, ribosomal activity, and many others) were found as part of the enriched routes for all the core genes, with functions that are in line with other similar pan-genome studies (Mosquera-Rendón et al., 2016;Valot et al., 2015). P. aeruginosa genome is composed of a mosaic structure including the large core-genome (Valot et al., 2015), into which regions of genomic plasticity lead to the insertion of a block of genes belonging to the accessory genome (Mathee et al., 2008). In the case of PaeAG1 and other ST-111 strains, the genome sequence is around 1.0 Mb longer that the reference genome Pae-PAO1, a difference that is re ected as genomic islands distributed along the genome.
Pae-RIVM-EMC2982 and Pae-Carb0163 (closest genomes to PaeAG1) had the most similar pro les carrying 41 out of the genomic islands. As highlighted in the Results, many genomic islands formed clusters (GICs, Fig. 3), including the genomic islands clusters harboring the two integrons (GIC VIM-2 and GIC IMP-18 ). Genomic island groups have been reported before as integrative and conjugative elements or ICEs (Petitjean et al., 2017), but ICEs in PaeAG1 (using ICEberg 2.0 platform, https://dbmml.sjtu.edu.cn/ICE nder/ICE nder.html) overlap with other GICs but none with GIC VIM-2 and GIC IMP-18 .
Since the size of the core-genome and its content is not well known (Valot et al., 2015), prediction methods are required to de ne accessory regions, but the outcome depends on algorithms (Ozer et al., 2014), which could explain differences and the GICs.
On the other hand, this prominent number of genomic islands in PaeAG1 and other ST-111 strains can be explained due to the absence of a functional CRISPR-Cas system (bacterial defense system against foreign DNA) and consequent high number of successful events of horizontal gene transfer (Petitjean et al., 2017). This genome plasticity of individual strains represents an advantage for P. aeruginosa to t the needs for survival in virtually any environment (Mathee et al., 2008).
In the context of carbapenems resistance, genes encoding for MBLs are usually found as gene cassettes in class 1 integrons (Jones-Dias et al., 2016;Walsh, 2005). This allows rapid dissemination in the clinical setting due to the selective pressure by the use of antibiotics (Sánchez-Martinez et al., 2010), which is aggravated since this antibiotic represents the last therapeutic resource to treat P. aeruginosa infections (Toval et al., 2015). While multiple studies correlate antibiotic resistance and the presence of integrons, the genetic context surrounding class 1 integrons is often not investigated in P. aeruginosa, as remarked before (Chowdhury et al., 2016).
Evaluation of the sequence showed that GIC VIM-2 is also present in Pae-RIVM-EMC2982 (100% coverage and 99.99% identity) and Pae-Carb0163 (100% coverage and 99.87% identity) at the chromosomal level.
However, a study including these strains showed that VIM-2-carrying integron and surrounding regions (~30 Kb, equivalent to GIC VIM-2 ) were shared with a plasmid of ST-446 P. aeruginosa S04-90 with 99% identity. Based on identity, mobilization of the fragment between plasmids and chromosomes may have occurred recently (van der Zee et al., 2018).
In the same study, analysis of genome landscape showed that the regions (equivalent to GIC VIM-2 ) corresponded to a DNA segment acting as a composite transposon, composed of four different transposons (Tn402, Tn21-like, a disrupted and another complete Tn4661). The class 1 integron carrying VIM-2 is contained in the Tn402 transposon (Gillings, 2017;van der Zee et al., 2018). Evolutionary details are completely explained in (van der Zee et al., 2018). GIC VIM-2 carries the genes involved in its transposition module (transposases such as TniB and TnpA) and mercury resistance module, as described in other similar transposons and insertion sequences (Chowdhury et al., 2016;Ghaly et al., 2017;Jones-Dias et al., 2016;Liebert, Hall, & Summers, 1999;van der Zee et al., 2018). The presence of gene cassettes unrelated to the antibiotic resistance can be a result of anthropogenic settings (Ghaly et al., 2017) and selection pressures in environments polluted with heavy metals and other substances such as mercury, arsenic, and disinfectants (Gillings et al., 2015).
Regarding the VIM-2-carrying integron, this element is an In59-like integron. In59 was rst reported two decades ago in France (Poirel et al., 2001) and then worldwide (Gillings, 2017;Samuelsen et al., 2010;Toval et al., 2015;van der Zee et al., 2018). Among all the 211 strains in our study, VIM-2 was only present in PaeAG1 and the two closest genomes (all ST-111). Differences in aacA29 genes de ned the aacA29e allele found in Pae-Carb0163 (van der Zee et al., 2018), in contrast to aacA29a-b in PaeAG1, all coding for aminoglycoside acetyltransferases. In the case of other VIM-2-carrying integrons, most of the 93 reported elements have an aacA allele, including the case of the In59-like present in PaeAG1 and other strains. A similar arrangement has been reported in other elements including In103, In927, In1025, and In1163 (van der Zee et al., 2018. Since GIC VIM-2 sequence and architecture are completely found in two VIM-2+/ST-111 strains, VIM-2carrying integron (In59-like) can be considered an old-acquaintance element in a well-known genomic context.
Additionally, the genomic context de ned by GIC IMP-18 was also analyzed. Using Pae-PAO1 as reference, it is shown that GIC IMP-18 insertion occurred in a speci c point (prrH) between PA4704 and PA4705 (Fig. 6).
This region contains three genes for regulatory small RNAs (prrF1, prrH, and prrF2) are found, which are involved in iron homeostasis under iron-depleted conditions (Reinhart et al., 2017) or to avoid iron toxicity (Reinhart et al., 2015).
While complete GIC IMP-18 (composed of GI48-GI49 genomic islands) was found in none of the other strains, the GI48 section was found in the Pae-97 strain (ST-234, with a class 1 integron), a genome close to the ST-111 group ( Fig. 2 and 3). Sequences comparison of GIC IMP-18 and Pae-97 showed 77% coverage and 99.92% identity. Gene composition of GIC IMP-18 includes endonucleases and recombinases module, the class1 integron, transposase TniB, and hypothetical proteins.
Since there is a lack of information about the genomic context of many IMP-carrying integrons (such as region GIC IMP-18 , unlike GIC VIM-2 ), and the particular architecture of the class 1 integron in PaeAG1 with the gene cassette IMP-18/gcuD/OXA-2/aacA4, we consider that this IMP-18-carrying integron (registered as In1666) is a novel element that we report here for the rst time.
In the partial reconstruction of the evolutionary steps related to the GIC IMP-18 region, the integron In1666, the insertion sequence IS1326, and the transposon TnAs3 seem to play a key role in the current state of this genomic region. Both IS1326 and TnAs3 have been reported in different integrons and high plasticity regions (He et al., 2016;Jones-Dias et al., 2016;Liebert et al., 1999;Szuplewska, Czarnecki, & Bartosik, 2014). Further analyses are required to complete the evolutionary steps which have de ned this genomic region as well as the implications of multiresistant in PaeAG1.
Jointly, identi cation of the landscape of the genomic context de ned by GIC VIM-2 and GIC IMP-18 , provides insights about the dissemination and evolution of mobile elements, in this particular case for integrons carrying MBLs. Since MBL-producing P. aeruginosa is associated with epidemic outbreaks and contributes to the dissemination of carbapenemase resistance worldwide (Castanheira, Deshpande, Costello, Davies, & Jones, 2014), it is worrisome that strains such as PaeAG1 circulate among Costa Rican hospitals. This can be correlated with the high prevalence of carbapenem-resistant strains in Costa Rica, many carrying VIM or IMP genes (Toval et al., 2015). Other ST-111 strains (harboring VIM-2 and IMP-18) from major Costa Rican hospitals are being characterized, and the genomic context for these isolates (including the possible presence of the new integron In1666) will take part in further studies to describe the possible evolutionary steps and the mobilization of the high-risk strains between hospitals. Future works are necessary to contribute the surveillance system to evaluate if other circulating strains carry these two elements, to identify its possible dissemination, and lead an adequate infection control program in medical centers.

Conclusions
PaeAG1 is a high-risk and critical organism due to its resistance to carbapenems by the activity of VIM-2 and IMP-18 enzymes, both harbored in two class 1 integrons. To describe the genomic context associated with these integrons, we rst veri ed the functionality of VIM-2 and IMP-18 after imipenem exposure. We then analyzed 211 complete genome sequences using a pan-genome analysis, separating strains by MLST pro le. Analysis of the 57 PaeAG1 genomic islands showed a varying pattern of the presence/absence among all the strains, in particular for closest genomes to PaeAG1. Two selected genomic island clusters, GIC VIM-2 and GIC IMP-18 , were studied in-depth. GIC VIM-2 sequence was completely found in the other two known ST-111 strains, which contained the VIM-2-carrying integron as an oldacquaintance In59-like element. GIC IMP-18 was partially found in another genome, but the IMP-18-carrying integron has an architecture never reported before, being considered as a novel In1666 integron. We provide new insights about the genomic determinants associated with this high-risk P. aeruginosa clone and its resistance to carbapenems using comparative genomics.

Availability of data and material
All the strains we used in this study were obtained from NCBI. All the IDs are available in the "Supplementary_ le 1 All_strains_information". For PaeAG1, the genome sequence and annotation les are available from our previous work at https://github.com/josemolina6/PaeAG1_genome, and in GenBank under the accession number CP045739.

Declaration of Competing Interest
The authors declare that there is no con ict of interest. * "Pa97_" is the locus with our annotation (see Methods). See Supplementary Table S2 for locus in PGDB annotation file and amino-acid comparison against other genomes. Cases with "CP913_" locus refers to the PGDB annotation file with a better score due to annotation algorithms differences. VIM-2 and IMP-18 expression after imipenem exposure. A RT-qPCR was performed to assess the transcriptomic activity of VIM-2 and IMP-18 genes. PaeAG1 was exposed to two imipenem concentrations, showing no effects on the growth curves (A). Relative gene expression showed that higher induction occurs at 6 hours after exposure, not only for VIM-1 (B) but also for IMP-18 (C). Relative expression was statistically different by time but not by concentration (p<0.05).

Figure 6
Description of the architecture of the exclusive genomic region GICIMP-18 containing the new IMP-18carrying integron. The genomic region GICIMP-18 is absent in the reference sequence Pae-PAO1 and Pae-RIV-EMC2982 strains, meanwhile it is partially present in Pae-97. The architecture of the IMP-18-carrying integron is shown with an arrangement that is reported here for the rst time.