Genomic insights into colistin resistance evolution in Pseudomonas Aeruginosa
- School of Biotechnology, International University, Ho Chi Minh City 70000, Vietnam
- Vietnam National University, Ho Chi Minh City 70000, Vietnam
- Research Center for Infectious Diseases, International University, Vietnam
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
() 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, 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 resistant to all tested antibiotics2. Colistin (polymyxin E) is an effective antimicrobial agent against gram-negative bacterial infections. However, it is often used as a last-resort treatment for infections caused by multidrug-resistant (MDR) and extensively drug-resistant (XDR) strains3. 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 charge6, 7. The other strategy is the overexpression of efflux pump systems8. Another mechanism is the overproduction of capsule polysaccharides9. 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 ATCC 9027 and analyze the genetic mutations caused by colistin exposure. By aligning the genomic sequences of colistin-exposed 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 underwent microevolution to survive antibiotic stress could be determined. With the emergence of colistin-resistant worldwide, a study on the development of colistin resistance in could be the key to controlling the spread of resistance.
Materials and Methods
Bacterial strains
ATCC 9027 was used as the original strain (PA-1). Strains Col-E1 and Col-E2 were obtained from ATCC 9027. In brief, 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.
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 data17. 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.
Workflow of the study.
General statistics of the
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.
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.
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 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 3A, B). The assembled N50 value was 284,316 bp, and the average GC content was 66.66% (Figure 3C). 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
The genomic alteration of 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 4A). A high overlap rate of variants between Col-E1 and Col-E2 was observed, including 153 SNPs and 448 INDELs (Figure 4B, 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
Functional annotation of Col-E1 and Col-E2.
|
Sample |
Type |
Impact |
Region Annotation |
Count |
Percent (%) |
|
Col-E1 |
SNPs |
LOW |
synonymous_variant |
45 |
1.67 |
|
MODERATE |
missense_variant |
27 |
1 | ||
|
HIGH |
stop_gained |
2 |
0.07 | ||
|
MODIFIER |
downstream_gene_variant |
133 |
4.94 | ||
|
MODIFIER |
intergenic_region |
88 |
3.27 | ||
|
MODIFIER |
upstream_gene_variant |
235 |
8.73 | ||
|
INDELs |
MODERATE |
conservative_inframe_deletion |
1 |
0.04 | |
|
MODERATE |
conservative_inframe_insertion |
3 |
0.11 | ||
|
MODERATE |
disruptive_inframe_insertion |
6 |
0.22 | ||
|
HIGH |
frameshift_variant |
14 |
0.52 | ||
|
HIGH |
frameshift_variant&start_lost |
1 |
0.04 | ||
|
MODIFIER |
downstream_gene_variant |
572 |
21.24 | ||
|
MODIFIER |
feature_elongation |
412 |
15.3 | ||
|
MODIFIER |
intergenic_region |
433 |
16.08 | ||
|
MODIFIER |
upstream_gene_variant |
721 |
26.77 | ||
|
Col-E2 |
LOW |
synonymous_variant |
45 |
1.66 | |
|
MODERATE |
missense_variant |
26 |
0.96 | ||
|
HIGH |
stop_gained |
2 |
0.07 | ||
|
MODIFIER |
downstream_gene_variant |
161 |
5.93 | ||
|
MODIFIER |
intergenic_region |
90 |
3.31 | ||
|
MODIFIER |
upstream_gene_variant |
264 |
9.72 | ||
|
INDELs |
MODERATE |
conservative_inframe_deletion |
1 |
0.04 | |
|
MODERATE |
conservative_inframe_insertion |
2 |
0.07 | ||
|
MODERATE |
disruptive_inframe_deletion |
1 |
0.04 | ||
|
MODERATE |
disruptive_inframe_insertion |
5 |
0.18 | ||
|
HIGH |
frameshift_variant |
13 |
0.48 | ||
|
MODIFIER |
downstream_gene_variant |
526 |
19.36 | ||
|
MODIFIER |
feature_elongation |
402 |
14.8 | ||
|
MODIFIER |
intergenic_region |
423 |
15.57 | ||
|
MODIFIER |
upstream_gene_variant |
756 |
27.82 |
A high number of alterations were found to have a MODIFIER effect (
As expected, the variant profile of Col-E2 was similar to that of Col-E1 (Figure 6,
Discussion
This study is the first in Vietnam to provide a whole-genome assembly of 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 environments19, 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 growth20, 21. Although the mechanism of polymyxin resistance is relatively well characterized22, there have not been many reports on genomic changes under sub-MIC antibiotic pressure. Thus, this article provides a preliminary overview of how the 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, , which results in the activation and overexpression of LPS-modifying genes 22. Among these, mutations in , and have been reported in colistin-resistant clinical isolates23.
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 in both Col-E1 and Col-E2. Mutations in are associated with resistance to streptomycin, kanamycin, and amikacin. Moreover, mutations were detected in many pathogens in the ESKAPE group24. Hence, under antibiotic pressure, the alteration of expression might be the key factor that helps pathogens adapt to extreme conditions. We also reported a mutation in the upstream region of , 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 was also found in in vitro-induced ciprofloxacin-, ceftazidime-, and meropenem-resistant strains26, 27. The genomic alteration in 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 28. Further investigation is needed to clarify the function of this mutation in the development of colistin resistance. In addition to , and are also important genes encoding the outer membrane proteins of . 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 LPS29. Nevertheless, with an deletion did not exhibit reduced susceptibility to polymyxin30. In contrast, the absence of OprF promoted biofilm formation and was linked to antibiotic resistance31.
Efflux pumps contribute significantly to multidrug resistance in 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 35. The efflux pump MexXY/OprM has been known to contribute to the resistance of to colistin36. In one study, the knockout of MexXY-OprM increased the susceptibility of , while mutations in and did not change the sensitivity to polymyxin37. Our study revealed an upstream mutation in which has transmembrane efflux transporter activity. This gene is known to help resist cobalt-zinc-cadmium38. 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 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 regions40, 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 and reported that these types of mutations increase or decrease the transcription of genes and are directly responsible for the evolution of important pathogenic phenotypes42. 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 levels43. 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 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 .
Conclusions
Overall, this study was the first in Vietnam to assemble the whole genome of 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 .
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.