Genomic insights into colistin resistance evolution in Pseudomonas Aeruginosa

Abstract

Colistin (polymyxin E) is the last-resort antibiotic used to treat multidrug-resistant Gram-negative bacterial infections, including those caused by Pseudomonas aeruginosa. The global emergence of colistin-resistant Pseudomonas aeruginosa is a significant concern for healthcare professionals. Therefore, identifying key contributors to the development of colistin resistance is crucial for addressing this issue. In this study, P. aeruginosa ATCC 9027 was serially exposed to subminimal inhibitory concentrations (MICs) of colistin for 14 consecutive days to obtain the Col-E1 strain. The Col-E1 strain was then cultured for 10 days in antibiotic-free medium to acquire the Col-E2 strain. The whole-genome sequences of three strains, namely, P. aeruginosa ATCC 9027 (the original strain, PA-1), Col-E1 and Col-E2, were assembled, annotated and analyzed. The bioinformatics pipeline included FASTQC (v0.11.9) for quality control; Unicycler (v.0.5.0) for de novo assembly; Bowtie2 (v2.4.5) and Picard Tools (v.2.27.4) for alignment; SAMtools (v1.11) and BCFtools (v.1.15) for variant calling; and SnpEff (v5.1) for variant annotation. As a result, we obtained a draft genome of P. aeruginosa ATCC 9027 consisting of 6,314,207 bp with 146-fold coverage. After aligning Col-E1 and Col-E2 against the draft genome, the number of insertion-deletions (INDELs) and single nucleotide polymorphisms (SNPs) were found, with 486 INDELs and 162 SNPs for Col-E1 and 474 INDELs and 163 SNPs for Col-E2. A high overlap rate of variants, including 448 INDELs and 153 SNPs, was observed between Col-E1 and Col-E2, indicating that the number of variants was constant. The analysis revealed notable mutations in genes encoding ribosomal components (rpsL, rpsG, and rpsJ), genes involved in efflux transmembrane transporter activity (czcC_1), and genes encoding an outer membrane porin (oprD). The variants were found in upstream (36.5%), downstream (25.7%) or intergenic (19.1%) regions. These mutations may be involved in gene expression regulation, leading to the development of a colistin-resistant phenotype. In conclusion, this study provided a preliminary overview of how P. aeruginosa responds to colistin antibiotic stress at the genomic level. These mutations could be used as markers for colistin resistance and further investigated to clarify their role in the colistin resistance of Pseudomonas aeruginosa.

Introduction

Pseudomonas aeruginosa ( P. aeruginosa ) is a gram-negative, rod-shaped bacterium that usually causes chronic infections. Due to its ability to adapt and survive under various environmental and physical conditions, P. aeruginosa is among the group of six bacteria known as “ESKAPE” pathogens and is considered a priority by the World Health Organization due to its antibiotic resistance 1 . In some cases, colistin is the only effective treatment for P. aeruginosa resistant to all tested antibiotics 2 . Colistin (polymyxin E) is an effective antimicrobial agent against gram-negative bacterial infections. However, it is often used as a last-resort treatment for P. aeruginosa infections caused by multidrug-resistant (MDR) and extensively drug-resistant (XDR) strains 3 . Colistin primarily acts on the outer membrane of gram-negative bacteria. It works by binding to the phosphate groups of membrane lipids and replacing divalent cations (Ca ²⁺ and Mg ²⁺ ). This process disrupts the integrity of the outer membrane and increases its permeability, causing the release of intracellular contents and ultimately leading to cell death 4 , 5 . Resistance to colistin occurs with lipopolysaccharide (LPS) modification through different routes. The most common strategies for overcoming colistin resistance are modifications of the bacterial outer membrane through the alteration of LPS and a reduction in its negative charge 6 , 7 . The other strategy is the overexpression of efflux pump systems 8 . Another mechanism is the overproduction of capsule polysaccharides 9 . In this study, whole-genome analysis using Illumina short-read sequencing and bioinformatic tools was carried out to assist in the characterization of antibiotic resistance. This study was the first in Vietnam to construct a full genome assembly of P. aeruginosa ATCC 9027 and analyze the genetic mutations caused by colistin exposure. By aligning the genomic sequences of colistin-exposed P. aeruginosa strains to those of the original strain before exposure, potential mutations that can be used as markers for colistin resistance could be identified, and how P. aeruginosa underwent microevolution to survive antibiotic stress could be determined. With the emergence of colistin-resistant P. aeruginosa worldwide, a study on the development of colistin resistance in P. aeruginosa could be the key to controlling the spread of resistance.

Materials and Methods

Bacterial strains

P. aeruginosa ATCC 9027 was used as the original strain (PA-1). Strains Col-E1 and Col-E2 were obtained from P. aeruginosa ATCC 9027. In brief, P. aeruginosa ATCC 9027 (colistin MIC= 4 µg/ml) was cultured in sub-MICs of colistin for 14 days to obtain Col-E1 (colistin MIC = 16 µg/ml). Col-E2 (colistin MIC = 8 µg/ml) was obtained after culturing Col-E1 in colistin-free medium for 10 days. The duration was selected on day 14 for exposure and 10 for reversion because they are the time points at which the MIC values of the strains became stable. All bacterial strains were stored at −80 °C in Tryptic Soya broth (TSB; HiMedia) supplemented with 30 % glycerol (TSB/glycerol 7: 3, v/v). For genomic DNA extraction, the samples were thawed directly from storage and cultured at 37 °C overnight at 120 rounds/minute (rpm) in TSB.

DNA isolation and whole-genome sequencing

Genomic DNA was extracted using a GeneJET Genomic DNA Purification Kit (Thermo Fisher Scientific). Total DNA was quantified using a BioTek Take3 (Agilent Technologies, Santa Clara, CA, USA). The paired-end (PE) libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, USA) according to the manufacturer’s instructions. Finally, the paired-end libraries (2 × 150 bp) were sequenced on an Illumina MiniSeq system at KTest Science Co., Ltd. (Vietnam). The raw sequences were deposited in the NCBI BioProject under accession number PRJNA1111315.

De novo assembly and sequence annotation

The PA-1 strain was used as a reference strain. The low-quality bases and adapter sequences were removed using Cutadapt (version 4.1) 10 . Sequence quality was evaluated using FastQC (v0.12.1) 11 . High-quality reads of PA-1 were then assembled by the hybrid assembly pipeline Unicycler (v.0.5.0), which yields accurate, comprehensive, and cost-effective data 12 . After that, the assembly was evaluated using QUAST (v5.2.0) 13 . Prokka (v1.14.5) was used to annotate the draft genome of PA-1 14 .

WGS sequencing analyses, processing, and variant calling and annotation

The raw sequence reads of Col-E1 and Col-E2 were also subjected to the QC process as described above. High-quality reads were mapped against the draft genome using Bowtie2 (v2.4.5) 15 with parameters specifying for paired-end reads followed by mark duplication using Picard tools (v.2.27.4). The alignment files were then sorted and indexed using SAMtools (v1.11) 16 . Qualimap (v.2.2.1) was used for evaluating sequence alignment data 17 . The SAMtools mpileup utility and BCFtools (v1.15) were used to call genomic variants. Hard filtering was then applied to the raw variant data to extract high-confidence mutations. High-confidence variants were subjected to SNPeff (v5.0e) 18 for gene-based and region-based annotation.

The workflow of the study is summarized in Figure 2 .

Figure 1 . Workflow of the study.

Figure 2 . General statistics of the P. aeruginosa draft genome. (A) Length in base pairs of assembled contigs. (B) Cumulative contig length. (C) GC content over contigs.

Figure 3 . Genomic variant profiles of Col-E1 and Col-E2. (A) Number of variants observed in Col-E1 and Col-E2. (B) Overlapping SNPs between Col-E1 and Col-E2. (B) Overlapping INDELs between Col-E1 and Col-E2.

Figure 4 . Functional annotation of the genomic variants of Col-E1. (A) Functional effects and (B) region-based annotation of SNPs. (C) Functional effects and (D) region-based annotation of INDELs.

Figure 5 . Functional annotation of the genomic variants of Col-E2. (A) Functional effects and (B) region-based annotation of SNPs. (C) Functional effects and (D) region-based annotation of INDELs.

Results

Draft genome sequence of Pseudomonas aeruginosa ATCC 9027.

After removing low-quality reads, 6,317,853 paired-end reads were retained, with an average quality per read of 36. A 6,314,207-bp assembly with 146-fold coverage was constructed consisting of 79 contigs ranging from 200 bases to 479,350 total bases in length ( Figure 3 A, B). The assembled N50 value was 284,316 bp, and the average GC content was 66.66% ( Figure 3 C). After draft genome annotation, 5980 genes were identified. The sequences ranged from 28 bp to 13,029 bp, including 5811 coding sequences, 99 miscellaneous RNAs, 69 tRNAs, and 1 CRISPR molecule.

Genomic variation of colistin-resistant Pseudomonas aeruginosa

The genomic alteration of P. aeruginosa in response to colistin was revealed by aligning both Col-E1 and Col-E2 against the PA-1 draft genome. High-quality mapping reads were obtained with 6,657,709 reads for Col-E1 and 7,717,242 reads for Col-E2. The mean depth coverages were 155.46-fold greater for Col-E1 and 181.52-fold greater for Col-E2, with a similar mapping rate of 99.38%.

After hard filtering, 162 and 486 high-confidence SNPs and INDELs were identified in Col-E1, while 163 and 474 high-confidence SNPs and INDELs were identified in Col-E2 ( Figure 4 A). A high overlap rate of variants between Col-E1 and Col-E2 was observed, including 153 SNPs and 448 INDELs ( Figure 4 B, C). This observation indicated that colistin exposure has a permanent effect on genomic alterations. Further investigations are necessary to reach a definitive conclusion.

Functional annotation of genomic variants.

The functional annotation of the genomic variants is summarized in Table 1 and visualized in Figure 5 and Figure 6 . In Col-E1, genomic alterations resulted in 2,693 functional effects, including 45 (1.67%) LOW, 37 (14%) MODERATE, 17 (0.63%) HIGH, and 2,594 (96.32%) MODIFIER. Among the HIGH effect variants, 2 stop_gain variants related to 2 SNPs at the ATCC_9027_59700_gene (p.Gln419*) and ATCC_9027_59770_gene (p.Arg65*) loci were identified. Most INDELs caused frameshift mutations 14 , and 4 INDELs were related to start_lost mutations ( Figure 5 , Table 1 ). These mutations are located in puuA and aroH, which are not known to cause antibiotic resistance.

Table 1 Functional annotation of Col-E1 and Col-E2.

A high number of alterations were found to have a MODIFIER effect ( Table 1 ). Notable genes included ribosomal components ( rpsL, rpsG, and rpsJ) , efflux transmembrane transporter activity ( czcC_1 ), and the outer membrane porin oprD . An insertion mutation upstream of oprD, which encodes an outer membrane porin that plays a significant role in the uptake of basic amino acids and carbapenems, was detected .

As expected, the variant profile of Col-E2 was similar to that of Col-E1 ( Figure 6 , Table 1 ).

Discussion

This study is the first in Vietnam to provide a whole-genome assembly of P. aeruginosa ATCC 9027 and characterize the genomic changes under colistin antibiotic stress. It has been suggested that antibiotics at low doses tend to increase the mutation rate in pathogens. The hypothesis is that higher mutation rates enable faster adaptation to specific environments 19 , 20 . Exposure to antibiotics has been shown to increase mutation and recombination frequencies in bacteria through the SOS response, even under sub-MIC conditions, which do not completely inhibit bacterial growth 20 , 21 . Although the mechanism of polymyxin resistance is relatively well characterized 22 , there have not been many reports on genomic changes under sub-MIC antibiotic pressure. Thus, this article provides a preliminary overview of how the P. aeruginosa genome responds to colistin antibiotic stress. The genomic profiles of Col-E1 and Col-E2 were similar, indicating that consistent genomic changes resulted in a resistance phenotype in both strains. However, the specific genomic variations that lead to the development of colistin resistance are still unclear.

As a last resort antibiotic, there are considerably fewer reports of colistin resistance compared to other antibiotics or polymyxins in general. The most common chromosome-encoded mechanisms of polymyxin involve mutations in two-component systems (TCSs), namely, PmrA/PmrB, PhoP/PhoQ, ParR/ParS, ColR/ColS and CprR/CprS , which results in the activation and overexpression of LPS-modifying genes 22 . Among these, mutations in pmrA , pmrB and phoQ have been reported in colistin-resistant clinical isolates 23 .

In this study, no common mutations were detected among polymyxin-resistant strains. This suggests that the development of colistin resistance involves regulatory mechanisms of gene expression. We identified insertion mutations occurring downstream of rpsL in both Col-E1 and Col-E2. Mutations in rpsL are associated with resistance to streptomycin, kanamycin, and amikacin. Moreover, rpsL mutations were detected in many pathogens in the ESKAPE group 24 . Hence, under antibiotic pressure, the alteration of rpsL expression might be the key factor that helps pathogens adapt to extreme conditions. We also reported a mutation in the upstream region of oprD , a porin that has been linked to many antibiotic-resistant strains, including carbapenem-resistant strains. Reduced expression/inactivation of OprD is a major carbapenem resistance mechanism 25 . Downregulation of oprD was also found in in vitro-induced ciprofloxacin-, ceftazidime-, and meropenem-resistant strains 26 , 27 . The genomic alteration in oprD found here might be linked to the alteration in the expression of OprD, which plays a part in colistin resistance in these strains. Indeed, the downregulation of OprD can also affect colistin resistance in P. aeruginosa 28 . Further investigation is needed to clarify the function of this oprD mutation in the development of colistin resistance. In addition to oprD , oprH and oprF are also important genes encoding the outer membrane proteins of P. aeruginosa . These mutations were detected, but no mutations were detected for these genes. OprH is thought to play an essential role in polymyxin resistance because it interacts with LPS and prevents polymyxin from binding to LPS 29 . Nevertheless, P. aeruginosa with an oprH deletion did not exhibit reduced susceptibility to polymyxin 30 . In contrast, the absence of OprF promoted biofilm formation and was linked to antibiotic resistance 31 .

Efflux pumps contribute significantly to multidrug resistance in P. aeruginosa 32 . Mutations in the efflux pump regions usually result in increased expression of the pump, which lowers the antibiotic concentrations inside the cell. These mutations are typically found in the promoter region or upstream of genes 33 , 34 . Mutations in the efflux pump MexA induce antibiotic resistance in P. aeruginosa 35 . The efflux pump MexXY/OprM has been known to contribute to the resistance of P. aeruginosa to colistin 36 . In one study, the knockout of MexXY-OprM increased the susceptibility of P. aeruginosa to cationic antimicrobial agents , while mutations in mexX and mexY did not change the sensitivity to polymyxin 37 . Our study revealed an upstream mutation in czcC_1, which has transmembrane efflux transporter activity. This gene is known to help P. aeruginosa resist cobalt-zinc-cadmium 38 . Although the connection between these heavy metals and colistin has been unclear, it has been reported that the combination of colistin with zinc oxide has a synergistic effect on P. aeruginosa 39 . The combination of colistin and metals can be further investigated, which could help increase the effectiveness of this last-resort antibiotic.

Similarly, the majority of variants appearing in the upstream and downstream regions suggest the role of regulatory factors in ensuring survival of the pathogen during the adaptation process. Previously, studies of pathogen adaptation and evolution have focused predominantly on coding regions 40 , 41 . However, noncoding regions (upstream, downstream, intergenic) might contain elements that regulate the expression of proteins and consequently determine virulence and resistance phenotypes. Khademi et al. investigated 44 clonal lineages of P. aeruginosa and reported that these types of mutations increase or decrease the transcription of genes and are directly responsible for the evolution of important pathogenic phenotypes 42 . More importantly, intergenic mutations help essential genes to become targets of evolution.

The development of antibiotic resistance is complex and involves changes at the genomic and proteomic levels 43 . These current characterizations of colistin-resistant genomes further suggest that the development of colistin resistance is a multifactorial process involving both changes to ensure survival and mutations that help individuals adapt to new environments. Further analysis of the resistomes and proteomes of colistin-resistant strains might provide a more comprehensive understanding of the development of colistin resistance.

This study provided an understanding of how P. aeruginosa with a fully antibiotic-susceptible profile and thus a basic genetic background respond to colistin exposure. Limitations of the study include the lack of repeated genomic sequencing of biological replicates to ascertain genetic changes due to colistin exposure and the lack of experimental evidence on the effect of each genetic change on phenotypic changes. Hence, further studies are required to clarify and improve our findings. Mutagenesis and mutational studies are essential to determine how these variants actually affect gene expression in the colistin resistance mechanism of P. aeruginosa .

Conclusions

Overall, this study was the first in Vietnam to assemble the whole genome of P. aeruginosa ATCC 9027 and characterize the genomic variants in in vitro-induced colistin-resistant strains. The analysis identified notable mutations that were mostly found in noncoding regions, indicating the possible gene expression regulation function of these mutations in the development of the colistin resistance phenotype. These variants could serve as potential colistin resistance markers for quick diagnostic tests. Mutagenesis and mutational studies should be further performed to understand the role of noncoding regions in the adaptation and evolution of colistin-resistant P. aeruginosa .

LIST OF ABBREVIATIONS

Whole genome sequencing (WGS), World Health Organization (WHO), multidrug-resistant (MDR), extensively drug-resistant (XDR), lipopolysaccharide (LPS), minimum inhibitory concentration (MIC), tryptic soya broth (TSB), paired-end (PE), single nucleotide polymorphism (SNP), insertion‒deletion (INDEL)

COMPETING INTERESTS

The authors declare that they have no competing interests.

ACKNOWLEDGEMENTS

This research was partially funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 108.04-2018.08.