Nucleotide substitutions in the mexR, nalC and nalD regulator genes of the MexAB-OprM efflux pump are maintained in Pseudomonas aeruginosa genetic lineages

Pseudomonas aeruginosa has different resistant mechanisms including the constitutive MexAB-OprM efflux pump. Single nucleotide polymorphisms (SNPs) in the mexR, nalC, and nalD repressors of this efflux pump can contribute to antimicrobial resistance; however, it is unknown whether these changes are mainly related to genetic lineages or environmental pressure. This study identifies SNPs in the mexR, nalC, and nalD genes in clinical and environmental isolates of P. aeruginosa (including high-risk clones). Ninety-one P. aeruginosa strains were classified according to their resistance to antibiotics, typified by multilocus sequencing, and mexR, nalC, and nalD genes sequenced for SNPs identification. The mexAB-oprM transcript expression was determined. The 96.7% of the strains were classified as multidrug resistant. Eight strains produced serine carbapenemases, and 11 strains metallo-β-lactamases. Twenty-three new STs and high-risk clones ST111 and ST233 were identified. SNPs in the mexR, nalC, and nalD genes revealed 27 different haplotypes (patterns). Sixty-two mutational changes were identified, 13 non-synonymous. Haplotype 1 was the most frequent (n = 40), and mainly identified in strains ST1725 (33/40), with 57.5% pan drug resistant strains, 36.5% extensive drug resistant and two strains exhibiting serin-carbapenemases. Haplotype 12 (n = 9) was identified in ST233 and phylogenetically related STs, with 100% of the strains exhibiting XDR and 90% producing metallo-β-lactamases. Haplotype 5 was highly associated with XDR and related to dead when compared to ST1725 and ST233 (RRR 23.34; p = 0.009 and RRR 32.01; p = 0.025). A significant relationship between the mexR-nalC-nalD haplotypes and phylogenetically related STs was observed, suggesting mutational changes in these repressors are highly maintained within genetic lineages. In addition, phylogenetically related STs showed similar resistant profiles; however, the resistance was (likely or partly) attributed to the MexAB-OprM efflux pump in 56% of the strains (only 45.05% showed mexA overtranscription), in the remaining strains the resistance could be attributed to carbapenemases or mechanisms including other pumps, since same SNPs in the repressor genes gave rise to different resistance profiles.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Haplotype 5 was highly associated with XDR and related to dead when compared to ST1725 and ST233 (RRR 23.34; p = 0.009 and RRR 32.01; p = 0.025). A significant relationship between the mexR-nalC-nalD haplotypes and phylogenetically related STs was observed, suggesting mutational changes in these repressors are highly maintained within genetic lineages. In addition, phylogenetically related STs showed similar resistant profiles; however, the resistance was (likely or partly) attributed to the MexAB-OprM efflux pump in 56% of the strains (only 45.05% showed mexA overtranscription), in the remaining strains the resistance could be attributed to carbapenemases or mechanisms including other pumps, since same SNPs in the repressor genes gave rise to different resistance profiles.

Background
Pseudomonas aeruginosa, a free-living microorganism, is a major opportunistic pathogen involved in healthcare-associated infections with high mortality rates [1,2]. Worldwide, P. aeruginosa is of major public health importance because it exhibits elevated resistance to nearly all types of antibiotics used for its mitigation [3,4] and the distribution of high-risk clones is increasing [5].
P. aeruginosa is known to have a non-clonal population structure with some outstanding sequence types (ST). The constant appearance of new variants has allowed for the observation that while this non-clonal structure is maintained in sensitive strains, multi-resistant strains primarily exhibit a clonal structure [5,6]. Worldwide, ST111, ST146, ST175, ST233, and ST235 are the major high-risk clones associated with multidrug resistance (MDR) and extensive drug resistance (XDR); other clones, including ST357 and ST664 have also been described but with less frequency [5,[7][8][9].
The MexAB-OprM efflux pump belongs to the resistance-nodulation-cell division family and is a complex that comprises three components: a cytoplasmic membrane transporter (MexB), a membrane fusion protein (MexA), and an external membrane pore (OprM). This system allows antimicrobials and other toxic compounds transportation from the cytoplasm to the extracellular environment [10,20].
Multidrug efflux pumps belonging to the RND family are mainly regulated at the transcriptional level, where regulators bind to DNA sequences that are located near or within the promoter regions of the genes they modulate [21]. Generally, transcriptional regulators contain one DNA-binding motif which consists of an α-helix-turn-α-helix motif. The second helix inserts into the major groove of target DNA and determine the binding specificity of the regulator; however, a third helix is needed to stabilize the DNA-binding motif [21]. In some cases, a local transcriptional regulator gene is linked to the same operon that encodes the efflux system it modulates, as occurs with mexR. In other cases, regulatory genes are dispersed throughout the chromosome and could regulate multiple genes simultaneously (nalC) [22].
Constitutive expression of the MexAB-OprM efflux pump maintains basal levels of the pump and is mediated by the repressor genes mexR, nalC, and nalD. The mexAB-oprM operon that encodes the efflux pump is regulated directly by the transcriptional repressors MexR and NalD and indirectly by NalC, which represses ArmR protein expression (MexR anti-repressor) to de-repress efflux pump expression. Mutations in the mexR, nalC, and nalD genes can impair their function, favoring MexAB-OprM efflux pump overexpression with a consequent increase in bacterial resistance [10,23,24].
Recent studies suggest SNPs in efflux pumps and porins are responsible for the generation of high-risk clones with MDR phenotypes [5,25,26]; however, in the case of the MexA-B-OprM efflux pump, only a small number of studies have investigated the relationship between SNPs in the three repressor genes (mexR, nalC, and nalD) and the MDR phenotype [27,28]. Mutational changes in the MexAB-OprM regulators may be associated with pressure exerted by the environment in which the strains are found, such as antibiotics, reactive oxygen species, detergents, among others [29,30]; however, it is unknown whether ST´s or high-risk clones of different origins present the same changes in these repressors regardless their genetic lineage.
This study aims to identify SNPs in the mexR, nalC, and nalD repressor genes of the MexA-B-OprM efflux pump in clinical and environmental isolates of P. aeruginosa (including highrisk clones) to discern whether these mutational changes may be associated with the pressure exerted by environmental conditions or are mainly related to genetic lineages.
Reference strains Escherichia coli ATCC25922, Pseudomonas aeruginosa ATCC 27853 (American Type Culture Collection, Manassas, VA, USA) and Klebsiella pneumoniae NCTC 13438 were used as controls to validate the different methodologies.
The strains were classified as follows: sensitive (S), resistant (R), multidrug resistant (MDR), extensively drug resistant (XDR), and pan drug resistant (PDR), according to the CLSI (2020) breakpoints and the criteria described by Magiorakos et al., 2012 [31, 32]. MDR strains were defined as non-susceptible to at least one agent in 3 or more antimicrobial categories; XDR, non-susceptible to at least one agent in all but two or fewer antimicrobial categories; PDR, non-susceptible to all agents in all antimicrobial categories.

Identification of carbapenemase-producing P. aeruginosa
Carbapenemase-producing P. aeruginosa were screened using the phenotypic technique βCARBA Test (BIO-RAD, France). This assessment is based on the color change of a pH indicator following hydrolysis of the β-lactam ring in carbapenem. For this test, all carbapenemresistant isolates were incubated at 37˚C for 24 h on Mueller-Hinton agar plates. Isolated colonies were recovered on a calibrated loop, resuspended in a bacterial protein extraction reagent, and further incubated for 30 min at 37˚C. Red color indicates a positive test, while no color change indicates negative. Escherichia coli ATCC 25922 was used as negative control, and Klebsiella pneumoniae NCTC 13438 (KPC-3 carbapenemase) served as the positive control [34].
For the mCIM test, a 10-μL loop of P. aeruginosa from an overnight blood agar plate was resuspended in 2 mL of Mueller Hinton broth, and a 10-μg meropenem disk was added following incubation for 4 h at 37˚C. A suspension of E. coli ATCC 25922 (0.5 McFarland, meropenem susceptible MIC � 2 μg/mL) was inoculated onto a Mueller-Hinton plate, and the previously-treated meropenem disk was added. For the eCIM test, 2 mL Mueller-Hinton broth with P. aeruginosa was prepared, and a 10-μg meropenem disk with 20 μL of 0.5 M EDTA was added and incubated for 4 h at 37˚C. Both meropenem disks were placed on one Mueller-Hinton plate inoculated with the meropenem susceptible reference strain and incubated for 24 h at 37˚C. The inhibition zones were measured and interpreted using the following guidelines established by CLSI, 2020 [32].

MexAB-OprM efflux pump phenotypic detection
Phenotypic activity of the MexAB-OprM efflux pump was assessed in the 91 P. aeruginosa strains and confirmed as a 4-fold decrease in MIC value for CB in the presence of Phe-Arg-βnaphthylamine (PaβN) inhibitor relative to that in the absence of the inhibitor. For this technique, the MexB-specific substrate CB was used as the reporter antibiotic [35]. The MIC was determined for each strain in the absence/presence of the efflux inhibitor PaβN (50 μg/mL). At this concentration, the inhibitor completely restored the susceptibility of the control strain P. aeruginosa ATCC 27853 (American Type Culture Collection, Manassas, VA) to the reporter antibiotic and did not inhibit bacterial growth [36,37]. All isolates grew in the presence of PaβN (50 μg/mL). The MexAB-OprM efflux pump was considered as the most likely cause of the elevation of the MIC (+) if the MIC value for CB -PaβN was at least 2 log 2 dilutions higher than in the reference strain (PAO1: MIC 64 μg/ml), and the MIC for CB +PaβN was lower than that measured in the reference strain (PAO1: MIC 64 μg/ml) or α 1 dilution; if the MIC values for CB +PaβN remained elevated compared with the reference strain, the MexA-B-OprM efflux pump was contributing in the elevation of the MIC ( � ); and if there was a difference of 1 dilution between CB -PaβN and CB +PaβN or no difference, the MexAB-OprM efflux pump was not the cause of the elevation of the MIC (-) [35].

MexAB-OprM efflux pump genotypic detection
P. aeruginosa from an overnight blood agar plate was resuspended in 2 mL of Mueller Hinton broth (0.5 McFarland) and incubated at 37˚C for 16-24h (synchronization phase) (each strain in duplicate). A 1:100 dilution was made, and bacteria was harvested at the late log-phase of growth. Total RNA was isolated using TRIzol Reagent (Zymo Research, USA) followed by DNase I treatment (Thermo Scientific, USA). Synthesis of cDNA was performed using the ZymoScript RT PreMix Kit (Zymo Research, USA) following the manufacturer´s protocol.
The  [35]. The efficiency of the reaction was calculated using the formula E = 10 (−1/slope) for which a standard curve was made with serial dilutions of cDNA (1:2) of the reference strain (PAO1) (in triplicate), the Cq values observed were plotted against the logarithm of the concentration to obtain the slope of each gene (mexA and rpsL). Gene transcript was normalized to the rpsL housekeeping gene and transcript levels were calculated as transcript expression ratio compared to the reference strain PAO1; the equation for relative quantification corrected for efficiency (Fold effect = (EmexA ΔCq,mexA(control−sample) )/ ErpsL ΔCq,rpsL(control−sample) ) was applied [27,39].

Bacterial DNA isolation and quantification
P. aeruginosa was cultured in Mueller Hinton broth at 37˚C for 18-20h. Chromosomal DNA was isolated using the Wizard Genomic DNA purification kit (Promega, USA), following the manufacturer´s protocol. DNA quality, integrity, and concentration were confirmed by agarose (1%) gel electrophoresis. DNA concentration and purity were evaluated using an EPOCH spectrophotometer (Biotek, Vermont, USA). High quality DNA was stored at 4˚C until used.

Genotyping via Multilocus Sequence Typing (MLST)
The ST of 58 of the P. aeruginosa strains was previously determined [6]; for the remaining 33 strains, the same genotyping procedure via MLST was performed. Nested PCR for the metabolic genes acsA, aroE, guaA, mutL, nuoD, ppsA, and trpE was carried out using the primers described by Curran et al., 2004 [40]. Sequencing of the PCR products was performed in both senses. The obtained sequences were edited and aligned as previously described. The ST of each strain was obtained by BLAST analysis (nucleotide) of each gene compared with the P. aeruginosa MLST database [41], http://pubmlst.org/paeruginosa/. The new STs were deposited in the P. aeruginosa MLST data base. Variability parameters were determined as previously described.

Phylogenetic analysis
The phylogenetic relationship and evolutionary relationship of the 91 nosocomial and environmental P. aeruginosa strains were evaluated by the construction of a phylogenetic network using maximum likelihood. To visualize the genomic relationships, a minimum-spanning tree was built from the MLST (ST) sequences using the GrapeTree [46], the phylogenetic inference was performed in PhyloViz Online [47], https://online.phyloviz.net/index softwares. In addition, groups of related STs (clonal complexes) were identified using the BURST analysis. A group of related STs was defined as a profile match at n-3 loci to any other member of the group (n = number of loci in the scheme, MLST = 7); default settings were used to achieve the most stringent definition. The GrapeTree, PhyloViz Online and BURST analysis softwares are available at the Pseudomonas aeruginosa PubMLST database [41], http://pubmlst.org/ paeruginosa/.
Furthermore, to determine the evolutionary relationships and events of recombination between the STs a phylogenetic network was built from the MLST (ST) of the P. aeruginosa strains using the neighbor-net algorithm (distance-based-method) implemented in SplitsTree ver.4.0. [48]. The robustness of the network was calculated with a bootstrap test after 1000 pseudo replicates and the inference of recombination events during the generation of allelic variation was estimated with the pairwise homoplasy index test (PHI).

Genetic diversity
For each of the MLST and MexAB-OprM efflux pump repressor genes, the number and frequency of haplotypes was determined, as well as the estimated nucleotide diversity, including the nucleotide diversity per site (Pi) and expected variation per site assuming a neutral evolution (Eta). The number of substitutions (S) for each gene is reported as well. All data were obtained using DnaSP 5.10.01 [45]. The DnaSP program allows the analysis of DNA polymorphisms using data from several loci by estimating several measures of DNA sequence variation within and between populations.

Statistical analysis
The nature of the data determined the type of statistical analysis used. Qualitative variables defined subgroups of the total cohort; therefore, associations between variables required the construction of contingency tables to identify association patterns from the counts within these values. Statistical significance was considered as p �0.05 as determined by Fisher´s exact test using STATA/MP 14.1 [49]. To investigate the effects of explanatory variables on a binary response variable we used logistic regression models, while for a categorical dependent variable with outcomes that have no natural ordering, multinomial logit models were used. All procedures were done with the STATA/MP 14.1 program [49].
For principal component analysis (PCA), a graphic was constructed using RStudio software [50], (http://www.rstudio.com/) to evaluate the relationships between variables (outcome, resistance, ST, MexAB-OprM haplotype, site, year, and ward). Through: (1) determining the covariance matrix of the normalized data, (2) finding the characteristic root and the characteristic vector, (3) determining the contribution rate of the variance of the principal components, (4) removing the main components, and (5) obtaining the principal component value and the integral score.
Statistical significance between mexA transcript expression and phenotypic detection of the MexAB-OprM efflux pump, susceptibility profile, and haplotype was evaluated by Kruskal-Wallis equality of populations rank test using the STATA/MP 14.1 program [49].

P. aeruginosa nosocomial strains exhibit more multidrug resistance than environmental strains
The susceptibility profiles of the studied strains to 14 antibiotics in nine categories is shown in Table 1. The 91 strains were classified as S, R, MDR, XDR, or PDR according to their susceptibility profile (Table 1) as determined by the CLSI 2020 standard values and criteria established by Magiorakos et al. [27,28]. 96.7% of the strains (88/91) were classified as multidrug resistant. Of the nosocomial strains, 49.35% (38/77) were classified as XDR; ten strains showed intermediate resistance to P/T, eight strains were AZT sensitive, and six strains were CS intermediate. In addition, 33.77% of the nosocomial strains were classified as PDR (26/77), 15.58% as MDR (12/77), and 1.30% (1/77) as sensitive. Of the environmental strains, 78.57% were classified as MDR (11/14), 7.14% as XDR (1/14), and 14.29% as sensitive (2/14). P. aeruginosa nosocomial strains were highly associated with the XDR profile compared to the environmental strains and taking the MDR strains as reference (RRR = 34.82; p = 0.001).
Carbapenemase expression was evaluated in 74 strains confirmed to be meropenem and/or imipenem resistant. The commercial kit β CARBA Test showed invalid results for 15 strains (Table 1). For the remaining 59 strains, carbapenemase typing (serine carbapenemase or    metallo-β-lactamase) was conducted ( Table 1) and showed that eight strains were positive for serine carbapenemases, and 11 strains were positive for metallo-β-lactamases.

Multilocus sequence typing verifies the diversity of most P. aeruginosa strains and reveals the emergence of outstanding sequence types
The ST for 58 P. aeruginosa strains were acquired from the PubMLST data base entry http:// pubmlst.org/paeruginosa/ (S1 Table), being the endemic clone ST1725 the most frequent and persistent for over 7 years in the hospital. The remaining 33 strains isolated during 2013-2015 were analyzed here, revealing three new alleles: allele 233 for aroE, 147 for guaA, and 157 for mutL; in addition to 23 new ST that were integrated into the worldwide P. aeruginosa MLST database (S1 Table). During the same period, six ST233 strains and one ST111 strain (both reported worldwide as high-risk clones) were isolated, with ST111 identified in the environment. Of the 48 ST studied, the environmental strains showed the greatest diversity, with a different ST in each strain ( Table 1). The nucleotide and gene diversity were greatest among environmental strains (Pi, 0.0074; Hd, 0.0071), with the greatest diversity observed in the aroE and trpE genes (Pi, 0.010; Hd, 0.007). Of the 94 SNPs identified, 81 were in hospital strains, with the greatest number seen in the aroE (n = 18) and trpE (n = 19) genes. Relevant genetic data, including the number of haplotypes, nucleotide diversity, gene diversity, and substitutions identified by MLST genotyping are summarized in Table 2.

SNPs in the regulatory mexR, nalC, and nalD genes show nalC gene as the most diverse
It was observed a total of 62 nucleotide substitutions, 49 synonymous and 13 non-synonymous in the mexR, nalC, and nalD genes (S2 Table) Haplotype was defined as the DNA sequence of the concatenated mexR-nalC-nalD genes. A total of 27 different haplotypes were identified in the 91 strains, including 26 haplotypes with substitutions and one haplotype with a 12-bp deletion (S2 Table and Table 3). The hospital strains showed the largest number of haplotypes (n = 19) (Table 1), while the environmental strains had the greatest diversity in nucleotides (Pi) and genes (Hd) (Pi, 0.00922; Hd, 0.879), with the greatest diversity observed for the nalC gene (Pi, 0.01184; Hd, 0.771) ( Table 2). We observed a total of 62 SNPs, of which 61 were in hospital strains and most were in the nalC gene (n = 28) (S2 Table).

Phylogenetic analysis arranges P. aeruginosa strains into genetic complexes that share the same characteristics including mexR-nalC-nalD haplotypes
The phylogenetic network based on the MLST genotyping (ST) of the P. aeruginosa strains is shown in Figs 1 and 2. Although both figures refer to the same phylogenetic network, different information is highlighted. The phylogenetic relationships, evolutionary relationship, clonal complexes identified in the 48 STs and the 27 mexR-nalC-nalD haplotypes are shown in Fig 1; in Fig 2, the relationship between the 27 identified mexR-nalC-nalD haplotypes, the 48 STs, the susceptibility profiles, and production of carbapenemases are depicted. The data indicated six important clonal complexes (CC); in each complex, all STs have a close phylogenetic relationship as follows: Complex I. Includes strains ST1725 as the most prevalent (n = 34), followed by ST1723, ST1730, ST2243, ST2244, ST2245, ST2247 and ST2248 (n = 1 strain each), all nosocomial origin (Fig 1). The ST1725 was identified by the BURST analysis as the potential Ancestral Type (AT) of this clonal complex (Fig 1). Globally reported data in the PubMLST P. aeruginosa database describe this complex as part of the CC309.
Haplotype 1 was highly associated with Complex I compared to other haplotypes and taking the singletons STs as reference (RRR = 409.53; p = 0.000). Haplotype 1 was identified with high frequency (n = 40); all ST1725 strains presented this haplotype except one isolate classified as haplotype 2 because of a deletion (4 105-116 ). ST111, a high-risk clone, also presented this haplotype (although it appears to be phylogenetically distant from CC1). Haplotype 1 was identified in strains of different ST, suggesting that these substitutions are specific to phylogenetically related ST (Fig 1).
Complex II. Includes ST1724, ST1726, ST1728 and ST1727, all nosocomial origin with the ST1724 identified as the AT (Fig 1). This complex form part of the globally CC235. All STs that make up this complex presented haplotype 5 and showed XDR (Figs 1 and 2). The relationship between this haplotype and fatal patient outcomes was remarkable. Complex II was associated with death when compared to the singletons STs (RRR = 40.01; p = 0.006), to complex I (RRR = 23.34; p = 0.009), and to complex 4 (RRR = 32.01; p = 0.025). The high-drug resistance exhibited by these strains could be in part attributed to the MexAB-OprM efflux pump activity; in two strains (2/5), efflux pump activity was the most likely cause of the elevation of the MIC, and in three strains (3/5) efflux pump was contributing to the elevation of the MIC. mexAB-oprM overtranscription was detected in two strains (Table 1 and Fig 2). Complex III. Includes ST1729, ST540, ST2250, and ST2251, all environmental origin except for one XDR ST1729 strain that was associated with a fatal outcome in 2009 (Figs 1 and 2). No AT was identified within these STs; however, they are part of the global CC253. All STs that make up this complex presented haplotype 26 and showed variations in antimicrobial susceptibility (S, MDR, and XDR) (Figs 1 and 2). In 2/6 strains, efflux pump activity was contributing to the elevation of the MIC. mexAB-oprM overtranscription was detected in two strains (1 MDR, 1 XDR) ( Table 1).
Complex IV. Includes ST2559, ST2560, and ST233 (n = 6), all nosocomial origin (Fig 1). ST233 is considered the AT of this complex, conforming the CC233 worldwide. All STs that make up this complex presented haplotype 12, were XDR and produced metallo-β-lactamases (Figs 1 and 2). In three strains (3/9), the activity of the pump was contributing to the elevation of the MIC. mexAB-oprM overtranscription was detected in seven strains (Table 1).
Complex V. Includes ST1737 and ST561, both of nosocomial origin and being part of the CC245 worldwide. Both STs that make up this complex presented haplotype 21 and are possibly associated with fatal patient outcomes (Fig 1). This complex presented MDR and the Mex-AB-OprM efflux pump was not the cause of the elevation of the MIC. mexAB-oprM overtranscription was not detected (Table 1 and Fig 2).
Complex VI. Includes ST2710, ST2704, ST2713, ST2716, and ST2731. All these STs were of nosocomial origin with ST2710 identified as the AT but no global CC identified. All STs that make up this complex presented haplotype 8, except for a ST2313 strain that presented  https://doi.org/10.1371/journal.pone.0266742.t003 haplotype 9, being the only difference between these haplotypes a SNP in the nalC gene S 209 R (Fig 1 and Table 3). This complex showed MDR (66.66% XDR and 33.33% PDR) (Fig 2). Only in one strain, the MexAB-OprM efflux pump was contributing to resistance; in the remaining strains, the pump was not the cause of the elevation of the MIC. Haplotype 8 strains produced serine carbapenemases and the haplotype 9 strain produced metallo-β-lactamases. mexAB-oprM overtranscription was detected in four strains (2 XDR, 2 PDR) ( Table 1). The rest of the STs are considered singletons since more than 3 locus variants are noticed between them and other STs. The remaining 19 mexR-nalC-nalD haplotypes were identified in the STs considered singletons (Fig 1). Most of these STs were MDR and only one strain (ST2568) produced serine carbapenemases (Table 1 and Fig 2).
Close phylogenetic relationship between CCI and CCII is evident, while the CCIV is the most distant complex. However, all CCs appear to be related somehow to the CCI (Fig 1).
In addition, the six clonal complexes previously identified in the phylogenetic networks stand out as groups of highly maintained STs in the neighbor-net graph (Fig 3). The presence of rectangular boxes in the network represents the high probability of extensive homologous recombination, which was corroborated with the PHI test that revealed statistically significant recombination events (p<0.05). According to the neighbor-net graph, CCII (Haplotype 5) and CCIII (Haplotype 26) appear to be closely associated with the CC1 (Haplotype 1), but distant from CCIV (Haplotype 12) and CCVI (Haplotype 8) which are close and are related to production of carbapenemases, indicating that despite the highly recombinant nature of P. aeruginosa, some substitutions are highly maintained among the strains (CCs). This close relationship between specific ST that differs little in sequence (CCs) is troubling, as it raises the possibility of the eventual selection of an efficient high-risk clone with high dissemination capacity.

Principal component analysis and statistical analysis reveal association between the ST and the mexR-nalC-nalD haplotype
PCA analysis showed a strong relationship between resistance and ST (Fig 4). The first two principal components of the analysis explained 66.15% of the observed variation, specifically the first component explained 47.81% and the second 18.34% (Fig 4). PCA also revealed In addition, the MexAB-OprM haplotype (mexR-nalC-nalD) were associated with the ST (p< 0.0001) and resistance (p< 0.0001) (this last one due to the association between ST and resistance).
PCA showed a relationship between the strain isolation site, isolation date, and hospital ward. However, these variables were inversely proportional to ST and haplotype. On the other hand, the ST and haplotype variables showed closeness and the same direction. The outcome variable (fatal outcome) was inversely proportional to all analyzed variables, although an association was observed between haplotype 5 strains (mexR-nalC-nalD) (Fig 4, green dots), haplotype 1 strains (blue dots), and fatal outcomes. Fatal outcomes were observed in 12.5% of haplotype 1 strains, 80% of haplotype 5 strains, 11.11% of haplotype 12 strains, and 16.67% of haplotype 26 strains (p = 0.051) (Figs 1 and 4, see haplotype colors). However, it should note that patients' underlying conditions were not considered in this study.
Statistical analysis of the relationship between the SNPs in the mexR-nalC-nalD genes, the sequence type, and patient death outcome gave the following results: 1. Sequence type (ST): It was observed that phylogenetically related sequence types presented equal or similar mexR-nalC-nalD haplotypes (p <0.05) (Fig 1, see Table). However, it should note that patients' underlying conditions were not considered in this study.

Discussion
This study identifies SNPs in the regulatory mexR, nalC, and nalD genes of the MexAB-OprM efflux pump in clinical and environmental isolates of P. aeruginosa (including high-risk clones) to discern whether these changes are associated with the pressure exerted by environmental conditions or are mainly related to genetic lineages. P. aeruginosa has a great capacity to resist adverse environments, as evidenced by its nosocomial survival. Worldwide, different studies have reported between 27.6 and 71.4% resistant isolates [6,24,51]. We observed a high fraction of PDR (33.7%) and XDR (78%) in strains isolated from a hospital setting, while those from environmental settings were MDR (78%). In our work and worldwide, the worrisome observation has been made that colistin resistant strains are present in both nosocomial and other environments while colistin is the last therapeutic option for MDR and XDR isolates [52,53]. Dößelmann et al. (2017) warn about the rapid acquisition of colistin resistance after 10 and 20 days of exposure [53]. The presence of antibiotics or heavy metals in the environment induces MDR and even XDR environmental strains. The finding of resistant bacteria in the environment is attributed to the discharge of antibiotics into waste waters [54] and in industrial waste and the misuse of antibiotics as a preventive measure in livestock and fish farms [4]. Regarding heavy metals, contamination could be accomplished to destruction of built; waste disposal; soil and rock erosion; and industrial and agricultural practices. Sub-inhibitory levels of heavy metals can produce selective pressure on bacterial populations contributing to MDR through co-selection mechanisms (co-resistance, cross-resistance, and co-regulatory resistance) [55-58]. These factors may explain our observed high resistance to numerous antibiotics (including fosfomycin) in environmental and clinical strains, which may interact in nature.
We identified a total of 48 different ST of P. aeruginosa strains of nosocomial and environmental origin (S1 Table). As expected, all the environmental strains had a different ST, representing great diversity and indicating a non-clonal population structure. In contrast, the nosocomial strains had a non-clonal population structure; however, the emergence of highly successful epidemic clones was evident [53] (ST1725 and ST233), as previously described by Aguilar-Rodea et al. [6].
Among the nosocomial strains, ST1725 stands out for its high frequency and multidrug resistance (PDR, 21 strains; XDR, 12 strains; MDR, 1 strain). This endemic clone (Complex I) prevailed for over seven years (2007)(2008)(2009)(2010)(2011)(2012)(2013) in the institution. Similarly, we identified six P. aeruginosa nosocomial strains as ST233 (Complex 4), all were XDR. Although ST233 has been previously reported worldwide as a high-risk clone resistant to AZT, there are no records of colistin resistant strains or even XDR or PDR strains [59-63]. ST233 is a risk for patients in Mexico, as these strains were resistant to colistin and only sensitive to AZT. The AZT susceptibility of our analyzed strains may result from the lack of commercial use of this antibiotic in Mexico due to restrictions. In the United States (Northeast, Ohio), where AZT use is free of such restrictions, AZT-resistant ST233 strains have been reported [60].
We also identified ST111 (Clomplex I) in a MDR environmental strain collected from a water source. ST111 is considered a high-risk clone identified in Croatia and France with MDR association [64-66]. Occasionally, P. aeruginosa strains migrate from the environment and cause animal or human infections [26,67]. For this reason, this type of ST is considered highly dangerous and should be kept under surveillance.
MDR P. aeruginosa strains likely result from several factors. The involvement of overexpression of efflux pumps has gained recognition, particularly the MexAB-OprM pump for its constitutive expression and attribution of resistance to most antibiotics [12][13][14][15][16]. Our results corroborate the relationship between the MDR of a given P. aeruginosa strain and its MexA-B-OprM efflux pump activity, principally in PDR strains (73%). Efflux pump activity was also higher in strains ST1725 and ST233, compared with other STs (p �0.0001).  [69]. In our study, of the P. aeruginosa strains that showed MexAB-OprM activity (phenotypic positive-strains, 56.04%), only in 52.9% the efflux pump was the most likely cause of resistance, in the remaining 47% of the strains the efflux pump was contributing to the resistance; p overtranscription was detected in 45.05% of the strains; the correlation between both methodologies was 52.75%; however, mexA transcript was detected in all the analyzed strains, verifying the MexAB-OprM efflux pump basal level expression. Similarly, Goli et al. (2018) reported a correlation of 66.6% and 68% between the results of MICs with PAβN and mexB gene overtranscription in clinical strains [19]. Other authors report that the high correlation between phenotypic and genotypic methodologies was only detected in reference strains or highly characterized laboratory strains, not in clinical isolates where diverse phenotypes were shown due to diverse genetic backgrounds [70][71][72]. Highlighting the difficult interpretation of phenotypic data due to the co-expression of resistance mechanisms other than efflux (enzymes, plasmid-encoded b-lactamases, mutations, class 1 integrons), and by the selection pressure exerted by treatment with various antibiotics on the patient; arguing that genotypic methods do not necessarily provide data on the final expression of the gene product and its functionality [19,70,71]. Molecular mechanisms such as post-transcriptional control of the protein translation rate, half-lives of proteins or mRNAs, RNA stabilization mechanism, and the molecular association of the protein products of expressed genes must be considerate [71,73]. Additionally, these regulatory mechanisms of transcription, translation, and proteolysis are highly affected by the genetic heterogeneity in clinical isolates, which can lead to discrepancies [71]. In this study, both methodologies demonstrated that although the MexAB-OprM pump shows positive activity and contributes to resistance, this is not the only mechanism of resistance in these strains [18].
SNPs in MexAB-OprM regulators, which are mainly attributed to environmental conditions, can impair their function favoring mexAB-oprM overtranscription expression [10,23,24]. Recent studies have identified point mutations in the mexR, nalC, and nalD repressor genes of the MexAB-OprM efflux pump [24,27]; however, it is unknown whether STs or highrisk clones of different origins present the same substitutions in these repressors regardless their genetic lineage. Here, we observed 27 different haplotypes in the three repressor genes. Regarding the mexR gene three non-synonymous substitutions were identified in 71.42% of the strains, with the V 126 E amino acid variation being the most frequent, agreeing with previous studies [27]. The 268 C!T nonsense substitution was the only change that encoded a stop codon (Q 90 � ) and was observed in two strains (haplotype 18). Finally, the 170 T!C nonsense substitution produced the amino acid variation L 57 P in the strain identified as 73H. Nevertheless, when the preliminary modeling of MexR was performed, apparently no effect was noticed on the structure of the protein when considering the amino acid variations V 126 E and Q 90 � ; however, a considerable conformational change was observed when modeling the amino acid variation L 57 P, but further studies are needed to know its effect on the final structure of the protein. Most of the non-synonymous substitutions were found in the nalC gene, as seen in 98.90% of the strains (90/91). The most frequent substitution was G 71 E (96.70%; 88/91) as previously described by Horna et al (2018) [27]. In addition, a 12-bp deletion was identified in one strain (S2 Table and Table 3). Another study also identifies the nalC gene as the main site of nucleotide substitutions target, reporting relevant changes (non-synonymous substitutions) in 87% of nosocomial isolates as well as some deletions [27]. Most of the substitutions reported in that study differ from those we observed (Table 3). Regarding the nalD gene, the only nonsynonymous substitution identified (A 46 T) was previously reported by Suresh et al (2018) [24]. A crystallography study of the NalD protein demonstrated 66.5% similarity and 32.5% identity to the TetR family protein (TtgR), with different binding ligands with antibiotics such as novobiocin (N129 and H167), in addition, the authors also reported that mutations in F175 and N176 alanine interferes in the binding of NalD to its promoter [74]; however, further studies are needed to know the effect of these changes on the final structure of NalC and NalD proteins, standing out that the crystallized structure for NalC has not been reported.
It has been reported that genetic variations as stop codons, frameshifts or deletions lead to loss of functionality of the repressor genes and may contribute to the mexAB-oprM over transcription [75] with a consequent increase in bacterial resistance [10,23,24]; however, in our study, same haplotypes (SNPs) lead to different MexAB-OprM efflux pump behavior and resistance. The 12-bp deletion identified in the nalC gene, could be one of the causes of the PDR in the haplotype 2 strain due to the malfunction of this repressor (mexAB-oprM overtranscription). In the case of the stop codon identified in the mexR gene, mexAB-oprM overexpression was detected, in addition to the production of metallo-β-lactamases, both mechanisms could explain the XDR shown by these strains. However, in some cases correlation between mexA transcript level and haplotype was not evident, may be because in these, other newly discovered regulators as MdrR1, MdrR2, BrlR, CpxR, ArmR and PA3225 could be potentially involved in the transcription of the pump genes [76].
Specific haplotypes were identified within the different complexes, suggesting that these substitutions are highly maintained within phylogenetically related ST no matter what environment they came from. In addition, high-drug resistance exhibited by some of these complexes, could be mostly attributed to the MexAB-OprM activity, however different behavior of the pump (mexA transcript expression) was noticed within strains with the same haplotype (for example: haplotypes 8 and 12 which correlate closely with XDR and PDR profiles, mexAB-oprM overtranscription in the 80% and 77.8% of the strains (respectively) was identified), suggesting that the association between haplotypes and high-drug resistance also may be due to potential relationship between the ST and mexR-nalC-nalD haplotypes, as resistance is highly related to specific STs.
As haplotype 12 was identified in all ST233 (high risk clones), the sequences (mexR-nalC-nalD) detected in this work were compared with six complete genome sequences reported for ST233 strains around the world [36,77], the finding show a 100% of similarity in the SNPs present in the MexAB-OprM efflux pump regulators which are maintained in P. aeruginosa genetic lineages regardless the strains origin. In addition, for these strains, the resistance could also be explained by other mechanisms, as explained by Correa et al. [8], the presence of carbapenemases stands out in haplotype 12 strains, as also described by other authors in ST233 strains worldwide [5]. Interestingly, the carbapenemase types differ between strains of haplotype 8 (serine carbapenemases) and those of haplotype 12 (metallo-β-lactamases) [5].

Conclusions
Pseudomonas aeruginosa is a highly recombinant microorganism with a great population diversity; however, core genome of this specie is highly conserved, there are genes that despite having a high selection pressure are conserved within genetic lineages, such as the housekeeping genes. For its part, the constitutive MexAB-OprM efflux pump which is related to resistance to antibiotics and virulence of the strains is regulated by three repressors genes. In this work it was demonstrated that nucleotide substitutions in the mexR, nalC and nalD genes are highly maintained in phylogenetically related strains despite the selection pressure exerted by the indiscriminate use of antibiotics. A clear example is the ST233 reported in other parts of the world, in which same mutations present in its repressors could be observed regardless of the geographic location of isolation.

Nucleotide sequence accession numbers
The nucleotide sequences obtained in this study were deposited in the GenBank database under the following accession numbers: mexR gene sequences, MT188163-MT188173; nalC gene sequences, MT188174-MT188190; nalD gene sequences, MT188191-MT188204.
Accession numbers for the Pseudomonas aeruginosa isolates used in this work are available at the public database for molecular typing: PubMLST.org. (See S1 Table) and CDC & FDA Antibiotic Resistance (AR) Isolate Bank (See S3 Table).