Genome-wide association study of signature genetic alterations among pseudomonas aeruginosa cystic fibrosis isolates

Abstract

Pseudomonas aeruginosa (PA) is an opportunistic pathogen that causes diverse human infections including chronic airway infection in patients with cystic fibrosis (CF). Comparing the genomes of CF and non-CF PA isolates has great potential to identify the genetic basis of pathogenicity. To gain a deeper understanding of PA adaptation in CF airways, we performed a genome-wide association study (GWAS) on 1,001 PA genomes. Genetic variations identified among CF isolates were categorized into (i) alterations in protein-coding regions, either large- or small-scale, and (ii) polymorphic variation in intergenic regions. We introduced each CF-associated genetic alteration into the genome of PAO1, a prototype PA strain, and validated the outcomes experimentally. Loci readily mutated among CF isolates included genes encoding a probable sulfatase, a probable TonB-dependent receptor (PA2332~PA2336), L-cystine transporter (YecS, PA0313), and a probable transcriptional regulator (PA5438). A promoter region of a heme/hemoglobin uptake outer membrane receptor (PhuR, PA4710) was also different between the CF and non-CF isolate groups. Our analysis highlights ways in which the PA genome evolves to survive and persist within the context of chronic CF infection.

Introduction

Pseudomonas aeruginosa (PA) is an ubiquitous gram-negative bacterium that can cause disease in plants and animals [1]. PA infection can cause acute syndromes such as pneumonia and bloodstream infections, as well as chronic airway infections in patients with cystic fibrosis (CF). CF is a well-known genetic disorder caused by a mutated cystic fibrosis transmembrane conductance regulator (CFTR) protein. CFTR disruption alters the condition of the lung such that the increasingly dehydrated viscous mucus layer provides a favorable habitat for several pathogens, such as PA, which can infect opportunistically and persist over an extended period of time [2].

The ecological versatility of PA is thought to be associated with its relatively large genome containing numerous regulatory genes that confer an advantage in adapting during prolonged infections [3]. Within-host mutation of PA isolates from a single CF patient have been analyzed to better understand the underlying genetics of adaptation [4]. Expression studies have also been performed on PA grown in CF sputum and have identified altered expression of genes encoding amino acid biosynthesis and degradation, and quinolone signaling [5,6]. The lasR gene, which encodes an important quorum-sensing regulator, has frequently been detected among isolates from CF patients [7,8], however, the significance of this mutation is unclear. Mutations leading to increased antibiotic tolerance have also been identified in PA inhabiting the CF airway [9–11]. Furthermore, PA clonal lineages vary across different CF patients, depending on the sputum composition [12–15], and even within an individual patient over time [16,17]. Together, these findings suggest that PA is an adaptable organism that responds flexibly to changing environments.

Most large-scale genetic studies have focused principally on CF isolates. Therefore, in the absence of a comparator non-CF PA population, it is difficult to identify specific genetic variations that are associated with the CF environment. In this study, we compared 1,001 PA genomes sampled from both CF and non-CF clinical isolates. Genome-wide association studies (GWASs) have been useful in uncovering causal relationships between genetic variations and disease phenotypes in human populations [18], with over 4,000 human GWAS conducted to date [19]. In contrast to human genomes, bacterial genomes of the same species can vary in size, gene repertoire, and gene arrangement [20]. Because of this genomic variability, additional considerations must be made when performing bacterial GWAS [21–25]. In this study, we selected the genomes of PA isolates of known origin and performed a GWAS based on k-mer counting [22,23,26], a modified method that permits association mapping of genes and intergenic regions of the PAO1 reference genome. Results provided here expand our current understanding of how gene-level changes correlate with mechanisms of PA adaptation to the CF lung environment.

Results and discussion

P. aeruginosa population structure

Most studies aimed at better understanding chronic PA infection in CF patients have focused only on CF isolates [14,16]. It is therefore unclear whether genetic variations highlighted in these studies are important in the context of adaptation to the CF lung. Here, we assembled a collection of 2,167 PA genomes from the Pseudomonas Genome Database [27], including genomes from both CF and non-CF individuals. Genomes of unknown origin were removed and a phylogeny was constructed using 1,001 genomes of known origin (Fig 1A). In order to remove clonal isolates that may increase any potential confounding lineage effects in the subsequent GWAS, the dataset was trimmed to 636 genomes (S1 Fig, S1 Data), maintaining 99.8% diversity of original selection and the overall structure of the phylogenetic tree. Isolates from CF patients could be identified from all parts of the phylogeny, including the previously described groups A and B; and the long branching groups C1 and C2 [28]. Most CF isolates (409 among 422) belonged to group A.

[Figure omitted. See PDF.]

Fig 1. Visualization of significantly different regions and proteins with large-scale mutation.

(A) Phylogenetic tree was constructed with 1,001 genomes with known host disease information. Black and red leaves each indicate non-CF and CF isolates, and the numbers of CF and non-CF genomes for constructing the phylogenetic tree are placed inside brackets. Most of the CF isolates are located in group A, which includes the majority of genomes, and the smaller groups are designated group B, C1 and C2. (B) 31-mers that were significantly associated with CF and non-CF groups were aligned to the whole PAO1 genome and visualized by Phandango. Regions marked by asterisks are ribosomal DNA sequences and ‘hs’ marks the hotspot region where a large number of candidate 31-mers align. The higher the point representing a 31-mer, the more significant the lrt p-value. (C) Proteins with large-scale mutation from the PAO1 genome and Lrt p-values calculated for clustering performed under two conditions (S95L95: 95% similarity and 95% length coverage; S70L95: 70% similarity and 95% length coverage) are described in the first and second columns as per the color index beneath the heatmap. ‘hs’ indicates candidate proteins located within the hotspot region in (B).

https://doi.org/10.1371/journal.ppat.1009681.g001

Identification of CF-associated genetic elements

In total, 636 genomes were used to perform a GWAS using Pyseer [29]. These comprised 206 genomes from CF and 430 from non-CF. A k-mer based approach was implemented to include variability within intergenic regions. We selected k-mers with a length of 31 base pairs (31-mer) consistent with other studies that have adopted similar approaches and tested the effect of changes in k-mer length on GWAS sensitivity and false positive associations [22,23]. Significantly associated k-mers (p-value cutoff = 2.8e-08) were aligned to the PAO1 reference genome and visualized using Phandango [30] (Fig 1B).

GWAS analysis identified 41,685 k-mers associated with CF with a p-value < 2.8e-08. Of these, 489 de-novo assembled contigs mapped to 109 protein-coding genes, six untranslated regions and 29 intergenic regions in the PAO1 reference genome. Four most meaningful peaks marked with asterisks in Fig 1B are rDNA sequences. A greater number of CF isolates did not seem to have 16S and 23S rDNA sequences compared with the non-CF group. However, we verified that several of these CF isolates (AU2342, AU6462, AU9739, AU10409, AU16821 and AU25116) do have 16S rDNA by PCR-amplifying and sequencing the 16S rDNA regions. By extension, we believe that the lack of detection of 16S and 23S rDNA sequences from several CF isolates have arisen from an in silico error. Despite this, the number of genes annotated using Prokka [31] were not significantly different between genomes and they were therefore included in the input database.

GWAS elements associated with allelic variation in P. aeruginosa from CF patients

Focusing on the k-mers that were mapped to the PAO1 reference genome, we obtained amino acid sequences for each gene from all isolate genomes. Two different homology clustering thresholds were used: 95% similarity and 95% gene length coverage (S95L95) and 70% similarity and 95% gene length coverage (S70L95) (S2 and S3 Data). GWAS results were then filtered by allele frequency using the likelihood-ratio test (lrt) p-values in Pyseer for both conditions (Fig 1C). This provided a list of candidate genes with multiple associated alleles. We hereinafter refer to these genes with multiple co-localized mutations to contain large-scale mutations. Variation in MucA (PA0763) served as a positive control, as MucA mutation is frequently observed in CF isolates [32,33].

Large-scale alterations in protein-coding regions

Most candidate associated genes were identified within the prominent hs region (Fig 1C), with the most significant lrt p-value (p-value = 2.02e-22) associated with a probable TonB-dependent receptor (S2 Data). Differences in the strength of association were noted across amino acid residues for PA2335 (Fig 2A). Analysis of raw filter p-values seemed to suggest an opposing conclusion in relation to the importance of PA2335 alteration in PA adaptation to the CF airway. However, these raw p-values are not weighted to account for the population structure (Fig 2B and 2C). Mapping variations in PA2335 against the phylogeny demonstrated alignment of the CF phenotype with the presence of its homolog, except for region a (Fig 2C). This region corresponds to the group B isolates (S1 Fig) and suggests that these isolates possess alternative CF associations. This pattern was also seen for variants in the PA2333, PA2334, and PA2336 genes. Variation in PA2332 was associated with the CF phenotype across the entire population suggesting it may have a broad influence in CF adaptation (S2 Data). This is consistent with a previous study that identified lineage-specific markers for groups A and B, including genes PA2333~PA2336 which were over-represented in group A [28].

[Figure omitted. See PDF.]

Fig 2. PA2335 candidate with the lowest lrt p-value in the large-scale mutation group.

(A) Total amino acid residues of PA2335 are presented and each amino acid locus is represented by individual squares. The color of each square indicates the lrt p-value assigned to that amino acid residue based on Pyseer results and colors are assigned according to the index above the heatmap, which was generated based on the negative logs of the maximum and minimum lrt p-values from the Pyseer result of amino acid locus. ‘N’ and ‘C’ mark the N- and C-terminus of PA2335, respectively. (B) Pyseer result of the PA2335 cluster is shown. Different conditions used in clustering and the frequency of the PA2335 large alteration are described in ‘condition’ and ‘af’ columns. Population structure is not considered in the calculation of filter p-values, but it is for the calculation of lrt p-values. (C) The phylogenetic tree of 636 genomes and their associated disease status are shown. Red and black lines in the ‘id’ row indicate genomes isolated from CF and non-CF patients, respectively. Black lines in the second and third rows indicate that the homolog of PA2335 clustered with the PAO1 protein at given clustering thresholds (S95L95:S70L95). Red lines indicate that the homolog did not cluster with the PAO1 protein. Genomes in region a are included in group B indicated in S1 Fig.

https://doi.org/10.1371/journal.ppat.1009681.g002

Functional annotation of proteins with large-scale mutations was performed using BlastKOALA [34]. KEGG ontology (KO) [35], definition, and pathway information of the candidate associated genes were identified in the CF or non-CF isolates (Table 1, S3 Data). Almost all of the large-scale mutations were more significantly associated with PA from CF patients. Large scale mutations of MucA (PA0763) are well known in CF isolates to result in increased alginate production and thus contribute to the mucoid phenotype, promoting infection [32,33,36,37]. PA2020 is a negative regulator of the MexXY efflux pump, which expels antibiotics from the cell [38]. Increased MexXY expression level and antibiotic resistance are induced by disruption to the PA2020 regulator. Moreover, CF sputum contains mucins, free amino acids, lactate, and glucose as a potential energy source and PA preferentially uses organic and amino acids over glucose as an energy source [39]. For this reason, we speculate that some CF isolates have evolved with mutations in glucose catabolism genes (PA2321 and PA3193) to avoid spending unnecessary energy in catabolizing a less optimal carbon source (glucose). Additionally, mannitol has been used as a therapeutic treatment to increase mucociliary clearance in the airway and simultaneously induce tobramycin sensitivity in PA by generating a proton motive force [40]. Therefore, we postulate that exposure through the concurrent therapeutic use of mannitol and antibiotics in CF patients may be the background from which the mutations arose in the gene cassette (PA2338, PA2340, PA2341, and PA2342). These large-scale mutations are highlighted in bold in Table 1.

[Figure omitted. See PDF.]

Table 1. KEGG annotation of candidate proteins with large-scale mutation.

All large-scale mutations were annotated by Blastkoala. Proteins with KEGG information are listed in the ‘Query’ column, and KEGG ontology (KO), definition, and related pathway information are shown. Proteins discussed in the results are highlighted in bold.

https://doi.org/10.1371/journal.ppat.1009681.t001

Small-scale alterations in protein-coding regions

Large-scale mutations are more likely to cause loss-of-function than scattered substitutions in homologous sequence (small-scale mutation). However, adaptation to the CF lung may also involve substitutions that alter protein-coding regions leading to additive effects. To investigate this, we analyzed the top 10 amino acid residues of the small-scale mutations (Tables 2 and 3). For the genes containing the most significant CF-associated variation, the functional genomics was further analyzed with laboratory validation experiments.

[Figure omitted. See PDF.]

Table 2. Top 10 amino acid residues of the small-scale mutation group with 31-mers aligned to the PAO1 genome.

Top 10 amino acid residues according to Pyseer results with 31-mers aligned to the PAO1 genome are presented. The reference protein corresponding to each residue is presented in the ‘protein’ column, and the ‘locus’ column contains the locus of the mutation and amino acid residue of the reference protein where the mutation is detected. SLR sequences of PA5438 are highlighted in bold. Multiple types of mutation at each locus of the reference protein may be present and detailed information is provided in S3 Data.

https://doi.org/10.1371/journal.ppat.1009681.t002

[Figure omitted. See PDF.]

Table 3. Top 10 amino acid residues of the small-scale mutation group whose 31-mers did not align to the PAO1 genome.

Top 10 amino acid residues according to Pyseer results whose 31-mers did not align to the PAO1 genome are shown. Reference proteins from genomes other than PAO1 and the PAO1 homolog of this reference protein are presented in columns ‘protein’ and ‘homolog’. The ‘locus’ indicates the mutation and amino acid of the reference protein where the mutation is detected. Multiple types of mutation at each locus of the reference protein may be present and detailed information is provided in S3 Data. The reference protein name in the “protein” column is composed of the genome id, connected by the latter underscore sign, to gene number within that genome. All sequences of the reference proteins in Table 3 are provided in S4 Data.

https://doi.org/10.1371/journal.ppat.1009681.t003

Mutations in a putative transcriptional regulator (PA5438)

Identification and experimental validation.

Mutations at amino acid loci 272–274 (SLR) in a probable transcriptional regulator gene (PA5438) were strongly associated with CF isolates (Table 2 and Fig 3A). The amino acid sequence ‘SLR’ was deleted in the 2,419th gene of the CF isolate 18A_661, a PA5438 homolog (Fig 3B). The SLR deletion in the PA5438 protein exhibits a highly positive correlation with CF isolates overall (Fig 3C), and the PA5438 sequences (except for the SLR sequence region) were exactly identical to the reference PA5438 sequence. Domain search with SUPERFAMILY 2 [41] predicted these amino acids to be a part of a sugar isomerase (SIS) domain (Fig 3D). Therefore, we constructed an in-frame deletion mutant, PA5438ΔSLR, to determine whether these three amino acids affect the transcriptional regulator function or not.

[Figure omitted. See PDF.]

Fig 3. SLR deletion in PA5438 increases the expression level of aceA and glcB.

(A) Total amino acid residues of PA5438 are presented and each amino acid locus is represented by individual squares. The color of each square indicates the lrt p-value assigned to that amino acid residue based on Pyseer results and colors are assigned according to the index above the heatmap, which was based on the negative logs of the maximum and minimum lrt p-values from the Pyseer result of amino acid locus. ‘N’ and ‘C’ mark the N- and C-terminus of PA5438, respectively. Residues at locus 272–274 of PA5438 with significant lrt p-values are underlined in red. (B) Comparison of the reference PA5438 to a homolog protein from a CF-isolated genome (18A_661_02419) is shown. The red arrow indicates deletion of SLR residues in PA5438 and corresponds to the red-underlined locus observed in A. (C) Phylogenetic tree of 636 genomes and their associated disease status are shown. Red and black lines in the ‘id’ row indicate genomes isolated from CF and non-CF patients, respectively. Red lines in the second, third, and fourth rows (272nd, 273rd, and 274th) each correspond to the absence of S, L, R residues, respectively. One exception to this is the red line in row 274 (asterisk-arrow), which was detected as R substituted by C. (D) A domain search result of PA5438 by SUPERFAMILY 2 is shown and the position of the SLR sequence inside the SIS domain is highlighted in bold. (E) AceA (isocitrate lyase) and GlcB (malate synthase) are enzymes involved in the glyoxylate shunt pathway (highlighted in red). (F) RNAs of PAO1, the PA5438ΔSLR mutant, and the complementation strain were extracted at OD600nm ~1.0 in LB medium and relative expression levels of aceA and glcB were measured. **p<0.01.

https://doi.org/10.1371/journal.ppat.1009681.g003

PA5438 is a transcriptional repressor that directly binds to the promoter region of the aceA (isocitrate lyase) gene and has been shown to repress the expression of the glcB (malate synthase) gene during growth in a non-C2 carbon source [42]. These repressed genes encoding enzymes involved in the glyoxylate shunt (GS) pathway (Fig 3E). To determine whether the suppressive activity of PA5438 is lost in the PA5438ΔSLR mutant, gene expression levels of aceA and glcB were measured by quantitative real-time PCR (qRT-PCR). The qRT-PCR results show that the expression level increased 17-fold for aceA and 5-fold for glcB (Fig 3F). Based on these findings, we expect that the deleted SLR sequence in PA5438 is an important region for determining aceA and glcB expression levels.

Further investigation of the PA5438ΔSLR mutant was performed to identify phenotypes that may aid in the adaptation of the strain to the CF environment. First, we assessed the phenotype of the mutant PA5438ΔSLR in LB media as growth rates of diverse CF isolates in LB media closely approximate those of CF isolates grown in artificial sputum medium (ASM) and synthetic CF sputum medium (SCFM) [17]. Slower growth in LB was observed in the PA5438ΔSLR mutant compared to PAO1 and the complementation strain, especially during the exponential phase (Fig 4A). In a recent study, RccR of Pseudomonas fluorescens, a PA5438 homolog, was shown to control pyruvate metabolism as well as the GS pathway [43]. Interestingly, pyruvate dehydrogenase (aceE) and aceA were regulated differently in the presence of 2-keto-3-deoxy-6-phosphogluconate (KDPG), which is bound to RccR [43]. We measured the relative expression level of aceE at the exponential phase because pyruvate metabolism affects growth. The expression of aceE in the mutant was suppressed compared to wildtype PAO1 and the complementation strain (Fig 4B). AceE generates acetyl-CoA from pyruvate, which feeds into the TCA cycle (Fig 4C). Therefore, we considered if the addition of metabolites related to this enzyme reaction (either as a precursor or product) would achieve recovery of growth. A defect in mutant growth was still observed with the addition of sodium pyruvate and sodium citrate, but the addition of sodium acetate as a substitute of acetyl-CoA resulted in a complete recovery of growth (Fig 4D).

[Figure omitted. See PDF.]

Fig 4. SLR deletion in PA5438 affects antibiotic susceptibility in LB medium.

(A) Growth curves of PAO1, PA5438ΔSLR mutant, and the complementation strain in LB were observed over 12 hours. **p<0.01; ***p<0.001 (B) RNAs of PAO1, PA5438ΔSLR mutant, and the complementation strain were extracted at OD600nm ~1.0 in LB and relative expression levels of aceE were measured. **p<0.01; ***p<0.001 (C) AceE (pyruvate dehydrogenase) is an enzyme involved in the entrance of The Krebs Cycle (highlighted in red). (D) Growth curves of PAO1, PA5438ΔSLR mutant, and the complementation strain in LB supplemented with 10mM sodium acetate were observed over 12 hours. (E) Antibiotic susceptibility tests against ciprofloxacin (Cip) and 10 mM sodium acetate (SA) were performed. Initial OD600nm of PAO1, mutant, and the complementation strain were adjusted to 0.005, and OD600nm was measured after 20 hours of static incubation in LB supplemented with Cip. The concentration of antibiotic was 0.125 μg/ml. ***p<0.001.

https://doi.org/10.1371/journal.ppat.1009681.g004

Since antibiotic tolerance can be associated with slow-growing or non-dividing “persister” bacteria [44], we measured susceptibilities of PAO1, PA5438ΔSLR and the complementation strain to different classes of antibiotics (tobramycin and ciprofloxacin) commonly used to treat PA infection in CF patients. In static cultures with tobramycin or ciprofloxacin, the PA5438ΔSLR mutant reached a significantly higher OD600nm in ciprofloxacin compared to wildtype PAO1 and the complementation strain, but only slightly higher OD600nm in tobramycin (Figs 4E and S2A). Since growth in LB recovered with sodium acetate (Fig 4D), we performed antibiotic susceptibility tests with sodium acetate supplementation. A limited increase of tobramycin susceptibility was observed in PA5438ΔSLR (S2A Fig), whereas susceptibility to ciprofloxacin was fully restored to that of the wildtype PAO1 and complementation strain (Fig 4E). Therefore, we postulate that slower growth is involved in the decreased susceptibility of the PA5438ΔSLR mutant to the two antibiotics.

Evidently, the composition of nutrients in LB media is very different to that of the CF sputum. Therefore, a synthetic medium was prepared to mimic the nutritional makeup of CF sputum (SCFM) and used it to measure growth. Unlike growth in the LB medium, in SCFM, wildtype PAO1, PA5438ΔSLR and the complementation strain showed similar growth (S2B Fig). Since the decreased expression level of aceE seemed to be an important factor of growth in LB medium, we measured the expression levels of aceE in PA cultured in SCFM. The expression level of aceE in the PA5438ΔSLR mutant was similar to that of the wildtype PAO1 and complementation strain (S2C Fig). In contrast, aceA and glcB expression levels were significantly increased in the PA5438ΔSLR mutant compared to the wildtype and complementation strain grown in SCFM (S2D Fig).

PA produces pyocyanin, an important virulence factor derived from phenazine and forms blue-green pigments. Pyocyanin induces reactive oxygen species (ROS) generation by transferring electrons to oxygen and increasing neutrophil apoptosis as a way of disrupting the host immune system [45]. As such, the supernatant color of a PA culture can function as a proxy for the virulent nature of the bacterium. Interestingly, we observed the 8 hour culture supernatant of PA5438ΔSLR in SCFM to be greener than that of the PAO1 and the complementation strain (Fig 5A). When we confirmed the relative pyocyanin levels of this bacterial culture, the PA5438ΔSLR mutant in SCFM resulted in 6~7 times more pyocyanin than PAO1 and the complementation strain (Fig 5A). Pyocyanin stimulates iron removal from transferrin in bacteria incubated under low oxygen conditions [46], and the level of transferrin is increased in the bronchoalveolar fluid from CF patients, compared to that of healthy individuals [47]. Therefore, we measured biofilms in which oxygen is limited, using the iron-saturated holo-transferrin. Compared to PAO1 and the complementation strain, the PA5438ΔSLR mutant produced less biofilm in SCFM, as demonstrated by a decreased OD550nm absorbance value. However, in the presence of transferrin, the mutant produced a similar level of biofilm whereas the PAO1 and the complementation strain produced less biofilms than in SCFM (Fig 5B).

[Figure omitted. See PDF.]

Fig 5. SLR deletion in PA5438 affects virulence and biofilm production in SCFM.

(A) After bacterial culture of PAO1, PA5438ΔSLR mutant, and the complementation strain in SCFM for 8 hours with shaking, a pyocyanin assay was performed with culture supernatant, and relative pyocyanin level was calculated using PAO1 pyocyanin level. Bacterial cultures used in pyocyanin assay are shown. **p<0.01 (B) After 32 hours of static culture of PAO1, PA5438ΔSLR mutant and the complementation strain in SCFM and SCFM supplemented with holo-transferrin (TF), biofilm was measured spectrophotometrically. **p<0.01; ***p<0.001.

https://doi.org/10.1371/journal.ppat.1009681.g005

Potential implications in CF airway infections.

These findings have potential implications in CF airway infections. An altered fatty acid profile has been associated with the airways of CF patients. For example, higher concentrations of palmitic acid and oleic acid have been detected in CF airway samples [48]. Fatty acid degradation upregulates the glyoxylate shunt (GS) pathway [49], and PA isolated from the CF lungs can induce the expression of genes involved in fatty acid metabolism and the GS pathway [50]. Furthermore, mucin is a major energy source in the CF lung environment in addition to fatty acids, and the GS pathway is important for mucin degradation and consumption [51]. Meta-transcriptomic analysis of several CF sputum samples in another study also found that aceA was upregulated and genes associated with glucose transporters and glycolysis were significantly downregulated [52]. Based on these previous studies, we assumed that the activation of the GS pathway may provide an advantage for PA proliferation in CF lungs.

Previous observations of GS pathway upregulation [50,52] are probably due to the presence of CF isolates that contain the abnormal PA5438 gene. In the PA5438ΔSLR mutant, aceA and glcB expression levels were increased not only in LB but also in SCFM. To investigate if the end products of the upregulated expression of aceA and glcB are functional, we examined the sequences of AceA and GlcB from genomes that harbor ΔSLR. Both promoters and proteins of AceA and GlcB were highly conserved (Table A, Table B in S1 Table). This indicates that the SLR deletion is the sole contributor to the transcriptional regulation of these two genes.

Glyoxylate shunt pathway is necessary for surviving phagocytosis and is upregulated in response to ROS [53,54]. We performed a gentamicin protection assay using our strains and mouse bone-marrow-derived neutrophils to measure how well the PA5438 mutant endures ROS-enriched and nutrient-poor phagocytic environment. We did not observe a significant difference in the CFUs of the PA5438 mutant strain compared to the wild type and the complementation strains in our experiment (S2E Fig).

As CF isolates have similar growth rates in LB and SCFM [17], we chose LB as the medium in which to investigate phenotypic differences between the wildtype and mutant PA. However, we observed a marked difference between the growth of the PA5438ΔSLR mutant cultured in LB and SCFM. In addition, the bacterial culture head a greener colour when SCFM was used (Fig 5A) than when LB was used. Given that PA5438 is involved in energy metabolism [43], we hypothesized that the mutant and wildtype phenotypes may vary depending on the type of medium used (LB or SCFM). Therefore, virulence phenotypes of the mutant, such as the increased pyocyanin levels observed in SCFM, are more likely to occur in the CF environment than the decreased antibiotic susceptibility observed in LB cultures. We also examined the genes involved in pyocyanin biosynthesis in the genomes that carry the PA5438ΔSLR mutation. Whereas several genomes were found to carry mutations in several of these genes when compared with the PAO1 sequences, other genomes possess a highly conserved gene set (S3 Fig). Based on this finding, we anticipate that virulent phenotypes such as that of the PA5438ΔSLR mutant are likely to manifest in the CF lung environment.

Despite differences in the nutrient compositions of LB and SCFM, the mutant phenotypes observed in LB cultures provide evidence that PA5438 acts as a transcriptional regulator. RccR of Pseudomonas fluorescens, a PA5438 homolog, has been shown to possess an SIS domain, the binding affinity of which is regulated by 2-keto-3-deoxy-6-phosphogluconate (KDPG) [43]. This KDPG-dependent moderation of PA5438 binding affinity is also dependent on the gene with which PA5438 binds. In the absence of KDPG, the binding affinity of PA5438 to aceE is increased, whereas its binding affinity to aceA is decreased [43]. A domain search using the SUPERFAMILY 2 database detected the SLR amino acid sequence of PA5438 to be inside the SIS domain (Fig 3D). The qRT-PCR results showed a decreased aceE gene expression level and an increased aceA gene expression level in LB (Figs 3F and 4B). Based on these findings, we propose that the SLR amino acid sequence may be a key region to which KDPG binds. In contrast to LB culture, in SCFM culture the aceE gene expression by the PA5438ΔSLR mutant is similar to that of wildtype PAO1 and the complementation strain. This observation may be attributed to the preferential use of amino acids over glucose in SCFM, which may lead to a decreased intracellular level of KDPG in the exponential phase of PA cultured in SCFM, unlike in LB. Further investigations are needed to conclusively determine whether KDPG indeed binds to the SLR region of PA5438.

Furthermore, elastase and pyocyanin are well-characterized virulence factors regulated by quorum sensing. We performed an elastase assay to test whether quorum sensing is related to the increased pyocyanin level [20]. Similar to the pyocyanin assay results, which showed a significant increase in pyocyanin in the PA5438ΔSLR mutant, relative elastase activity was slightly increased in the mutant compared to wildtype PAO1 and the complementation strain (S2F Fig). Therefore, we anticipate the increased pyocyanin level to be mediated by the quorum-sensing system. Pyocyanin stimulates iron removal from transferrin when PA is grown in low oxygen conditions [46], and increased levels of transferrin have been observed in CF patients [47]. Biofilm is a matrix of extracellular polymeric substances which protects the bacterial cells beneath, and oxygen supply is limited in developed biofilms [55]. Biofilm production of the PA5438ΔSLR mutant is decreased compared to PAO1 and the complementation strain in SCFM. In SCFM with transferrin, the mutant maintains the level of biofilm production, whereas PAO1 and the complementation strain produced less biofilms (Fig 5B). This implies that the mutant maintenance of biofilm production level in the presence of transferrin may confer greater fitness in the transferrin-enriched CF environment. Based on our experimental results, we speculate that this mutant may be problematic in the CF lung due to its increased virulence and biofilms.

Mutations in the L-cystine transporter (YecS, PA0313)

Identification and experimental validation.

A major limitation of aligning 31-mers to the PAO1 reference genome is that the 31mers that do not align due to insertion mutations are overlooked. To determine the reference genes for such insertion mutations, un-aligned 31-mers were de novo assembled, and the contigs that formed were mapped to 635 genomes (excluding PAO1). The top 10 amino acid residues with corresponding 31-mers that do not align to the PAO1 genome are described in Table 3. Three amino acid residues at loci 162–164 (Fig 6) had significant lrt p-values compared to other amino acid residues. At these loci, the YecS homolog had insertions of three amino acid residues, SLI (Fig 6). Moreover, there was an additional copy of SLI in the AU17965_3981_04951 protein, comprising a total of three stretches of ‘SLI’, in contrast to a total of two stretches in the PAO1 YecS (S4 Fig). SLI insertions were strongly associated with CF isolates (Fig 6C) and the YecS homologs of these isolates were identical to the reference AU17965_3981 sequence, except for the additional SLI in PAO1 YecS. The AU17965_3981_04951 protein is a cytoplasmic-membrane transporter, and the additional ‘SLI’ sequence is thought to span both cytoplasmic and transmembrane regions (S5A and S5B Figs) [56].

[Figure omitted. See PDF.]

Fig 6. SLI insertion in YecS reduces L-cystine utilization activity.

(A) The amino acid residues of CF-derived AU17965_3981_04951 protein are presented and each amino acid locus is represented by individual squares. The color of each square indicates the lrt p-value assigned to that amino acid residue based on Pyseer results and colors were assigned according to the index above the heatmap. The color index was generated based on the negative logs of the maximum and minimum lrt p-values from the Pyseer result of the amino acid locus. ‘N’ and ‘C’ mark the N- and C-terminus of AU17965_3981_04951, respectively. Residues presenting especially low lrt p-values are marked with a red underline and were found at loci 162–164 of the AU17965_3981_04951. (B) Comparison of the reference AU17965_3981_04951 to its PAO1 homolog YecS (PA0313) is shown. The red arrow indicates insertion of SLI residues in AU17965_3981_04951 and corresponds to the red-underlined locus observed in (A). (C) Phylogenetic tree of 636 genomes and their associated disease status are shown. Red and black lines in the ‘id’ row indicate genomes isolated from CF and non-CF patients, respectively. Red lines in the second, third, and fourth rows indicate the presence of SLI insertion in PAO1 YecS. (D) Growth curves of PAO1, YecSΩSLI mutant and the complementation strain cultured in M9 minimal media with MgCl2 instead of MgSO4, supplemented with glucose and L-cystine as the sole sulfur source were measured over 30 hours. **p<0.01; ***p<0.001 (E) After bacterial culture of PAO1, YecSΩSLI mutant, the complementation strain and the ΔpchΔpvd mutant, a pyochelin- and pyoverdine- deficient strain, in SCFM for 6 hours with shaking, the pyochelin in the bacterial culture supernatant was measured. *p<0.05; **p<0.01.

https://doi.org/10.1371/journal.ppat.1009681.g006

To experimentally assess if the SLI insertion effects a phenotypic change, an in-frame mutant was constructed with the additional SLI sequence inserted into the corresponding locus of PAO1 YecS. No difference in growth was observed between the wildtype PAO1 and the YecSΩSLI mutant when both strains were grown in SCFM (S5C Fig). However, in M9 minimal media supplemented with glucose and L-cystine as the sole sulfur source, the mutant grew slower than PAO1 and the complementation strain, whereas growth in M9 media with an additional sulfur source present (MgSO4) showed no difference (Figs 6D and S5D).

Two L-cystine transporters are known in E. coli, an ATP-binding cassette (ABC) importer (FliY-YecSC) and symporter YdjN [57]. Previous work has shown that the presence of either of the two transporters is sufficient to meet the L-cystine requirement for normal growth when L-cystine is used as the sole sulfur source [57]. We performed a blastp search of 636 genomes using the database of 238 YdjN homologs [58] and concluded that 635 out of 636 genomes do not seem to encode a YdjN homolog (S5 Data). Therefore, based on our blast search and experimental result (Fig 6D), we hypothesize that PAO1 only possesses the ABC importer system for L-cystine uptake. Moreover, the decreased growth of the SLI insertion mutant is attributed to the decreased activity of the ABC transporter rather than a complete loss of function. A comparison of growth between a YecS clean deletion mutant and our SLI insertion mutant would help in better understanding the assortment of L-cystine transporters present in PA.

Potential implications in CF airway infections

In addition to the SLI insertion mutation, we focused on other patterns of mutation in YecS that were observed across the CF isolate genomes. One such pattern is the deletion of an SLI sequence that results in CF isolates encoding just one SLI sequence in YecS homologs (S4A Fig). Another pattern of mutation observed by multiple alignments of YecS homologs is the deletion of long stretches of amino acid residues (S4B Fig). Since the SLI insertion mutant exhibited decreased L-cystine transporter activity (Fig 6D), we anticipate that such large deletions in the YecS protein would incur a complete loss of function.

Amino acid concentrations in the sputum vary across individual CF patients [59] and cysteine levels can be low or even undetectable [59]. Since cysteine is oxidized to cystine, a low level of one may indicate a low level of the other in the CF lung. Therefore, on one hand, large deletions in YecS may provide a survival advantage by allowing the CF isolate to avoid the energy cost of maintaining a functional YecS transporter in the cysteine-deficient CF lung environment. On the other hand, the reduction of intracellular cystine to cysteine contributes to ROS production in E. coli [60]. Thus, we expect small-scale mutations in YecS may be a result of evolution that protects the CF isolates from ROS stress in the cysteine-present CF lung environment. Additionally, pyochelin is a well-known siderophore produced by PA. Pyochelin may inflict tissue damage in the CF lung by inducing a continuous inflammatory response [61]. Interestingly, cysteine is a precursor molecule required for the production of pyochelin [62]. Therefore, we tested the effect of SLI insertion in YecS on pyochelin production of PAO1, the YecSΩSLI mutant, the complementation strain, and a ΔpchΔpvd mutant [63] that does not produce pyochelin and pyoverdine. The SLI insertion mutant exhibited decreased pyochelin production compared to the wildtype and the complementation strain (Fig 6E). A decreased level of pyochelin may be associated with a less virulent phenotype of the mutant, which may be an adaptation strategy for establishing a chronic infection in the CF lung environment. Whereas several genomes carried mutations in pyochelin biosynthesis genes, when compared with the PAO1 sequences, other genomes possess a highly conserved gene set (S6 Fig). Therefore, it may be the case that less virulent phenotypes such as that of the YecSΩSLI mutant proliferate in the CF lung environment.

Small-scale alterations in intergenic regions: phuS and phuR

Comparison of genomes using k-mers enables the investigation of intergenic regions. Of all contigs assembled from significantly different 31-mers, between CF and non-CF isolates, 29 contigs aligned to the non-coding regions of the PAO1 genome. Most aligned to either the hs region or rDNA sequences (Fig 1B). Since mutations located in the hs region may cause gene deletions, and rDNA sequence regions were excluded from analysis due to difficulties in interpretation, we selected regions other than these loci. Among several such intergenic regions, one intergenic region between phuR and phuS operons, was involved in the pseudomonas heme utilization (phu) systems (Fig 7A). The function of such a system is the acquisition of iron from the heme group of hemoglobin [64]. To evaluate the potential role of these mutations in CF adaptation, Pyseer analysis was performed using the nucleotide sequences of this intergenic region from PAO1 and other isolates (Fig 7B). The most significant hit was the transition of the 117th residue on the forward strand from cytosine (C) to thymine (T) (Fig 7B). Two mutations (C117T and C122T) are included in the phuR promoter region (from -35 region to transcriptional initiation site (+1)) (Fig 7C).

[Figure omitted. See PDF.]

Fig 7. One SNP in the phuR promoter region increases expression level.

(A) The 180 bp intergenic region between phuR and phuS is shown. (B) The top 5 associated SNPs detected using Pyseer within this intergenic region are listed and the number in the ‘residue of mutation’ column indicates the forward strand-based locus of mutation. The nucleotide in front of the locus number indicates the reference base in the intergenic region, whereas the nucleotide following the locus number indicates the changed base. Information contained in the following three columns is similar to those described in Fig 2B. (C) Promoter region of phuR (region from -35 to +1) based on the PRODORIC database [84] is visualized with WebLogo [85], and SNPs listed in (B) are marked by asterisks. The height of a stack reflects the degree of sequence conservation at that position, while the height of symbols within a stack indicates the relative frequency of each nucleotide at that position. The upper sequence is derived from the intergenic regions of CF genomes and the sequence below is based on intergenic regions of non-CF genomes. (D) Red and black lines in the ‘id’ row indicate genomes of CF and non-CF isolates, respectively. Red lines in the bottom five rows indicate mutations at each locus, shown in (B), and yellow lines indicate that this intergenic region is missing from the genome. (E) Culture of phuR::lacZ and phuRC117T::lacZ grown to OD600nm ~0.25 were used in a β-galactosidase assay, to compare phuR promoter activity. ***p<0.001.

https://doi.org/10.1371/journal.ppat.1009681.g007

Consistent with our findings, the frequency of mutations within this intergenic region in a CF isolate were found to increase significantly compared to the background mutation rate, and the mutated intergenic region increased phuR promoter activity [65]. In some cases, the intergenic region was deleted (Fig 7D) and the promoter activity is probably absent in these isolates. When such cases of promoter deletion are excluded from the analysis, a stronger correlation is expected between the mutations of this region and CF isolates. To examine whether the C117T point mutation affects the expression of phuR, we performed a promoter activity assay by measuring β-galactosidase activity. Promoters with this mutation exhibited increased activity (Fig 7E). Therefore, increased phuR promoter activity will ultimately result in increased phu system activity. We speculate that the C117T mutation is an important strategy for iron uptake and thus survival in the CF environment.

Conclusions

In this study, we compared a large set of genomes from clinical PA isolates (CF vs. non-CF) and identified mutations that occurred differentially in either CF or non-CF isolates using GWAS based on 31-mer counting (Fig 8). To integrate these findings within a functional genomics context [66], the impact on bacterial growth- and virulence-associated phenotypes was also assessed. Importantly, PAO1-derived variants, PA5438ΔSLR and YecSΩSLI, exhibited distinct phenotypes in in vitro assays. Moreover, a single nucleotide replacement in the promoter region of the phuR gene caused a considerable increase in gene transcription. Together, these results suggest that PA may benefit from small-scale mutations when establishing chronic infections in the CF airway. Mutations in the lasR gene encoding a QS regulator have been reported as frequent in CF isolates [7,8]. However, we found no evidence that mutations in LasR or other QS-related proteins were significantly overrepresented in CF or non-CF isolates. Given that the MucA mutation was represented in our analysis, demonstrating the precision of our bioinformatic approach, the lasR mutation might in fact be a common feature in PA, either causing chronic CF infection or other acute infections. Finally, while our study improves understanding the genetic changes in PA associated with chronic CF airway infection, they do not reveal the order in which these genetic changes occurred. Evaluating the chronology of adaptative changes in the CF lung will be important for explaining the persistence of PA.

[Figure omitted. See PDF.]

Fig 8. Schematic diagram with candidate genes.

Gene products of both (1–4) large- and (5–7) small-scale mutations significantly found in CF are illustrated with BioRender.com. We anticipate several phenotypic changes to arise as result of these mutations: (1) Large-scale mutations in the MucA protein stimulate alginate production which is an important component of biofilm; (2) PA2321 and PA3193 are involved in glucose catabolism. Disruption of PA2321 and PA3193 may be results of evolution in the CF environment which is enriched in free amino acids, nutrients preferred by P. aeruginosa; (3) Large-scale mutations in the PA2020 regulator upregulates expression level of MexXY, which may in turn result in increased antibiotic resistance; (4) Mutations in mannitol utilization gene cassette may aid in defence against the combined treatment of aminoglycoside and mannitol; (5) SLR deletion in PA5438 induces the glyoxylate shunt pathway, pyocyanin production and maintenance of biofilm in the presence of transferrin; (6) SLI insertion in YecS decreases cystine uptake and pyochelin production; (7) Mutation in phuR promoter region increases phuR expression.

https://doi.org/10.1371/journal.ppat.1009681.g008

Experimental procedures

Bioinformatic analysis.

Overall, 2,187 PA genomes were downloaded from the Pseudomonas Genome Database [27], and 1,001 genomes with information on the host disease were selected. A phylogenetic tree was created with RapidNJ, using Snippy-generated alignments of the core genome as input [67,68], and visualized with Microreact [69]. Based on the initial tree, the selection of genomes was trimmed down with Treemmer [70] to include 636 isolates while maintaining 99.8% of the original diversity (S1 Data). The prediction of protein-coding genes in the 636 genomes was performed using Prokka [31]. K-mers of 31 base pairs length (31-mers) were counted in the 636 individual genomes with Fsm-lite [71]. A similarity matrix was constructed with Snp-sim [72] using a core alignment of 636 genomes generated by Snippy [68]. Pyseer [29], which uses a mixed model (FaST-LMM), was run using the counted k-mers and the similarity matrix as input, and lrt p-values were assigned to each 31-mer. To determine the 31-mers that were significantly different between CF and non-CF groups, 31-mers with poor chi values (indicative of invalid ??2 test) or lrt p-values above 2.8e-08 (cutoff value in this analysis) were removed. As a result, 41,685 31-mers were detected to have significant associations with either disease status (CF or non-CF). Distribution of these 31-mers aligned to the PAO1 whole genome was visualized with Phandango [30]. Subsequently, de novo assembly of these 31-mers was completed using Trinity [73]. Overall, 489 contigs were constructed, and these contigs mapped to 29 intergenic regions, six untranslated regions, and 109 translated genes of the PAO1 genome with Blastn. Individual contigs that aligned to either protein-coding regions or intergenic regions were blasted against 635 genomes (the PAO1 genome was excluded) [74]. The top hit from each pairwise alignment with a blast e-value <0.01 was chosen as a homolog of the candidate sequence in each isolate.

In 453 contigs mapping to protein-coding regions, homologs were clustered under several conditions of similarity and gene length coverage (S95L95 -: 95% similarity and 95% length coverage; S70L95–70% similarity and 95% length coverage) using Blastclust [74]. Homologs that clustered with the reference PAO1 protein were assigned as ‘1’ and those that did not were assigned ‘0’. We then used the clustering matrix as input for Pyseer. Clusters with either poor chi values or lrt p-values larger than 3.22E-05 were removed. In this manner, genes with multiple co-localized mutations were determined and defined as ‘large-scale’ mutations (Fig 1C). Other protein-coding region mutations were treated as ‘small-scale’ mutations.

Where contigs mapped to protein-coding regions or intergenic regions, multiple alignment of the candidate sequence and its homologs was performed using Mafft [75]. This allowed characterization of nucleotides (for intergenic region candidates) or amino acid sequences (for protein-coding genes) at each locus for detecting small-scale mutations. If the residue of a homolog was identical to that of the reference PAO1 intergenic region or protein, it was assigned a value of ‘1’ at that locus, and if not, a value of ‘0’ was assigned. In this way, a locus matrix was created for use as the input for Pyseer and loci with either poor chi values or lrt p-values larger than 4.01E-06 were removed.

Finally, k-mers that failed to align to the PAO1 genome were also analyzed. However, the PAO1 reference gene in these cases was replaced with the reference gene detected in alternative isolate genomes. As such, we were able to investigate the presence of insertion mutations. For KEGG functional annotation and pathway analysis, BlastKoala [34] was used with large-scale mutations. Sequence comparison was performed by Clustalw [76] and multiple alignment was visualized with Jalview [77]. SUPERFAMILY 2 [41] search was used to investigate the domains of the PA5438 gene. Transmembrane domain prediction of the AU17965_3981_04951 protein was carried out using Phobius [56].

Bacterial strains and growth conditions.

All bacterial strains and plasmids used in this study are shown in Table C in S1 Table. PAO1 was used as a reference strain [1]. Bacterial cultures were grown in LB medium (1% [w/v] tryptone, 0.5% [w/v] yeast extract, and 1% [w/v] sodium chloride) at 37°C and Synthetic Cystic Fibrosis Medium (SCFM) using a previously described protocol [59]. Single bacterial colonies were picked from LB plates and inoculated in fresh LB or SCFM broths for precultures and grown overnight. Precultures were diluted 100-fold in fresh LB broth to subculture and incubated at 37°C with shaking at 230 rpm. For SCFM cultures, bacterial preculture was diluted and inoculated to the starting OD600nm 0.025 in fresh SCFM broth for subculture. All SCFM cultures were incubated under the same conditions as LB cultures, except for the incubation time, which was dependent on specific experimental procedures.

In-frame mutant construction

All in-frame mutants, including an SLI insertion in the YecS protein and SLR deletion in the PA5438 protein, were constructed from the PAO1 strain. In-frame insertions and deletions were performed to construct mutants with amino acid level changes. In the case of the in-frame deletions, the 5′ and 3′ flanking regions were designed to overlap in both directions of SLR. For the in-frame insertion, nucleotide sequences corresponding to amino acid sequence ‘SLI’ were inserted into the middle of the 5′ flanking region. An overlap was constructed by using each 5′ primer of 5′ flanking region and 3′ primer of 3′ flanking region. For both deletion and insertion mutations, an overlapping product was inserted into modified pCVD442, a suicide vector, containing gentamicin and ampicillin resistance markers. PAO1, grown on LB agar, was conjugated with E. coli SM10 λpir, harboring pCVD442 with the overlapping product inserted, grown on LB agar with 50 μg/ml ampicillin and 30 μg/ml gentamicin. Conjugates were spread onto LB agar with 50 μg/ml gentamicin and 20 μg/ml irgasan to select single crossover recombinants. The single crossover recombinant was incubated on LB agar without NaCl but containing 8% sucrose for the selection of the desired mutant. Sequence verification was performed by PCR. Primers used in constructing these mutants are listed in Table D in S1 Table.

Complementation

Identical primers were used in constructing the aforementioned mutants (Table D in S1 Table; PA5438#1, PA5438#4, yecS#1, and yecS#6), and were used to amplify the yecS and PA5438 sequences of wildtype PAO1. The amplified wildtype sequences were each inserted into pCVD442, and identical steps of mutant construction were followed to generate the complementation strains of the YecSΩSLI and PA5438ΔSLR mutants.

Promoter assay

To perform the phuR gene promoter assay, intergenic regions (180 bp) between phuS and phuR in the PAO1 genome and AU2342_3932 [78] with only the 117th locus changed from cytosine to thymine were amplified with primers listed in Table D in S1 Table. These PCR-products were each cloned into the upstream region of the β-galactosidase gene of puc18-mini-Tn7t-Gm-LacZ [79] for chromosomal insertion. The constructed plasmid was transformed into E. coli DH5α λpir. After verification by DNA sequencing, the plasmids with the helper plasmid pTNS2 that encodes TnsABC+D genes (allowing Tn7 transposition) were electroporated into PAO1. Empty puc18-mini-Tn7t-Gm-LacZ plasmid with no insert was used as a control to measure the baseline expression of lacZ. The potential clones were selected on LB agar with 50 μg/ml gentamicin and sequence verification was performed to select the final candidates whose transposon was inserted properly into the region following the glmS gene. β -galactosidase activities of the three clones (con::lacZ, phuR::lacZ, and phuRC117T::lacZ), summarized in Table C in S1 Table, were grown in LB and measured at the exponential phase (OD600nm ~0.25).

Growth curves

Wildtype PAO1, the PA5438ΔSLR mutant, and the complementation strain were precultured overnight in LB broth. Precultures were diluted 100-fold in fresh LB and incubated at 37°C with shaking at 230 rpm. Growth in LB was observed over 12 hours and the OD600nm was measured. Growth in SCFM was detected in the same manner over a period of 8 hours. Growth recovery tests of the PA5438ΔSLR mutant were performed using LB medium and LB supplemented with 10 mM sodium acetate (Sigma-Aldrich, USA), 10 mM sodium pyruvate (Sigma-Aldrich, USA), and 10 mM sodium citrate (Sigma-Aldrich, USA), and OD600nm were measured at regular time intervals. To examine the functional capacity of the L-cystine transporter, overnight precultures of PAO1, the YecSΩSLI mutant, and the complementation strain were washed in phosphate-buffered saline (PBS), and diluted to OD600nm ~ 0.025 in fresh M9 minimal media supplemented with 22.2 mM glucose and L-cystine (0.25 mM) as the sole sulfur source (2 mM MgCl2 instead of MgSO4). OD600nm was measured over 30 hours. Growth was also measured over 21 hours under identical conditions, except for the inclusion of 2 mM MgSO4 (instead of MgCl2) in M9 minimal media.

Antibiotic susceptibility test

For antibiotic susceptibility tests, ciprofloxacin (Duchefa, The Netherlands) and tobramycin (Sigma-Aldrich, USA) were used. Wildtype PAO1, the PA5438ΔSLR mutant, and the complementation strain were precultured overnight in LB broth. Bacterial preculture was diluted in LB containing the relevant antibiotics (ciprofloxacin or tobramycin) and the starting OD600nm was adjusted to 0.005. For static cultures, the antibiotic concentrations used were 0.125 μg/ml ciprofloxacin and 1 μg/ml tobramycin. OD600nm was measured after 20 and 26 hours of static incubation with ciprofloxacin and tobramycin, respectively. An antibiotic susceptibility test with 10 mM sodium acetate supplementation was conducted under the same conditions.

Reverse transcription and quantitative real-time PCR

PAO1, the PA5438ΔSLR mutant, and the complementation strain were precultured and subcultured in LB. After incubating the subcultures to OD600nm ~1.0, RNeasy Mini kit (Qiagen, Netherland) and on-column DNase1 digestion were used to extract RNA following the manufacturer’s protocol. 1 μg of RNA was reverse-transcribed to synthesize complementary DNA by using reverse transcriptase (Takara Bio, Japan) and random hexamer primers. To check for DNA contamination, the same cDNA synthesis was performed in the absence of reverse transcriptase. SYBR green-based qPCR was performed using an ABI 48-well StepOne real-time system. The primers used are listed in Table D in S1 Table. CT values were normalized by 16S rRNA CT values. SCFM cultures (OD600nm ~0.9) followed identical RNA extraction and qRT-PCR procedures.

Elastase and pyocyanin tests

Precultures of PAO1, the PA5438ΔSLR mutant, and the complementation strain were prepared for elastase and pyocyanin tests in SCFM as described above. Bacterial supernatant was harvested from bacterial subcultures incubated for 8 hours in SCFM, and assays were performed as described previously [80,81]. Elastase and pyocyanin values were normalized using OD600nm.

Pyochelin measurement

Precultures of PAO1, the YecSΩSLI mutant, and the complementation strain in SCFM were prepared as described above. Bacterial supernatant was harvested from bacterial subcultures, incubated for 6 hours in SCFM, and fluorescence was measured using Varioskan Flash 3001 (Thermo Scientific, USA) at excitation wavelength 355 nm and emission wavelength 440 nm [82]. Pyochelin values were normalized using OD600nm.

Biofilm assay

Precultures of PAO1, PA5438ΔSLR mutant, and the complementation strain in SCFM were prepared as described above. Bacterial precultures were diluted in SCFM or SCFM containing the holo-transferrin (Sigma-Aldrich, USA) in a 96-well plate, and the starting OD600nm was adjusted to 0.025. For static cultures, the holo-transferrin concentration used was 1.25 μM. After 32 hours of static incubation, biofilm assay was performed as described previously [83].

In vitro gentamicin protection assay

Neutrophils were isolated from the bone-marrows of C57BL/6 mice following the protocol of EasySep Mouse Neutrophil Enrichment Kit (STEMCELL, Canada). Following isolation, neutrophils were maintained in Opti-MEM media (ThermoFisher Scientific, USA). Precultures of PAO1, the PA5438ΔSLR mutant, and the complementation strain in SCFM were prepared in the manner described above. Bacterial culture was harvested from bacterial subcultures incubated for 4 hours in SCFM. The initial infection dosage was 107 CFU per 5 ✕ 105 neutrophils. After 2 hours of co-culturing P. aeruginosa and neutrophils, the culture supernatant was removed and fresh Opti-MEM supplemented with 50 μg/ml gentamicin was added for 1 hour. After performing wash twice in PBS, 0.5% Triton-X was added and bacteria number were counted by using pseudomonas isolation agar (Sigma-Aldrich, USA).

Statistical analysis

Data are expressed as mean ± standard deviation. Unpaired Student’s t-test (one-tailed, unequal variance) was performed to analyze the differences between experimental groups. P-values smaller than 0.05 were considered statistically significant. All experiments were repeated for reproducibility.

Supporting information

S1 Fig. Phylogenetic tree constructed with CF and non-CF genomes.

Upper phylogenetic tree was constructed with 1,001 genomes containing host disease information, and the tree below was drawn using 636 genomes and maintaining 99.8% diversity of the upper tree. Black and red leaves each indicate non-CF and CF isolates, and the numbers of CF and non-CF genomes for constructing each phylogenetic tree are placed inside brackets.

https://doi.org/10.1371/journal.ppat.1009681.s001

(TIF)

S2 Fig. Phenotypes of the PA5438ΔSLR mutant.

(A) Antibiotic susceptibility test with tobramycin (Tob) and Tob with 10 mM sodium acetate (SA) were performed. Initial OD600nm of PAO1, PA5438ΔSLR mutant, and the complementation strain were adjusted to 0.005, and OD600nm was measured after 26 hours of static incubation in LB supplemented with Tob. The concentration of antibiotic was 1 μg/ml. **p<0.01 (B) Growth curves of PAO1 and PA5438ΔSLR mutant in SCFM were observed over 8 hours. (C) RNAs of PAO1, the PA5438ΔSLR mutant, and the complementation strain were extracted at OD600nm ~0.9 in SCFM and relative expression levels of aceE were measured. (D) Relative expression levels of aceA and glcB of the same RNA used in (C) were measured. ***p<0.001 (E) The initial infection dosage was 107 CFU per 5 ✕ 105 neutrophils isolated from the bone-marrows of C57BL/6 mice. After 2 hours of co-culturing P. aeruginosa and neutrophils, we added gentamicin to the culture medium for 1 hour to remove any extracellular bacteria. We then harvested the intracellular bacteria by treating the neutrophils with 0.5% Triton-X, and measured the bacterial CFU. (F) After bacterial culture of PAO1, PA5438ΔSLR mutant, and the complementation strain in SCFM for 8 hours with shaking, an elastase assay was performed with culture supernatant. **p<0.01.

https://doi.org/10.1371/journal.ppat.1009681.s002

(TIF)

S3 Fig. Blastp results of genes associated with pyocyanin biosynthesis.

Each row represents a genome containing an SLR deletion in PA5438 homologs and color indicates whether it is included in CF or non-CF. Columns contain PAO1 proteins associated with pyocyanin biosynthesis and similarity (S), % of identical matches between reference PAO1 protein and its homolog, and length coverage (L), % of reference PAO1 protein sequence covered by its homolog. Gray box indicates there is no homolog under e-value 0.01 in the blastp search.

https://doi.org/10.1371/journal.ppat.1009681.s003

(TIF)

S4 Fig. Multiple alignment of YecS and its homologs.

Multiple alignment of YecS and its homologs is visually represented. Names of the genomes are shown to the left of the multiple alignment. PAO1 and AU17965_3981 are representative genomes of non-CF and CF groups, and genomes labeled with red or blue belong to the CF or non-CF groups, respectively. Black arrow marks where the additional SLI insertion occurs (162nd to 164th residues) compared to the YecS protein. Regions marked by asterisks are regions of SLI amino acid repeat sequences in YecS. Genomes in a contain a deletion of SLI, resulting in a single copy of SLI, and genomes in b present large deletions in YecS.

https://doi.org/10.1371/journal.ppat.1009681.s004

(TIF)

S5 Fig. Transmembrane domain prediction and growth curves of the YecS mutant.

(A) Predicted transmembrane domains of the AU17965_3981_04951 protein are portrayed. Numbers below the figure indicate the amino acid loci. (B) Detailed amino acid ranges of the predicted transmembrane domains are listed. The SLI insertion in the AU17965_3981_4951 protein is present within the region highlighted in bold. (C) Growth curves of PAO1 and YecSΩSLI mutant in SCFM were measured over 8 hours. (D) Growth curves of PAO1 and YecSΩSLI in M9 minimal media supplemented with glucose, L-cystine, and MgSO4 were recorded over 21 hours.

https://doi.org/10.1371/journal.ppat.1009681.s005

(TIF)

S6 Fig. Blastp results of genes associated with pyochelin biosynthesis.

Each row represents a genome containing an SLI insertion in YecS homologs and color indicates whether it is included in CF or non-CF. Columns contain PAO1 proteins associated with pyochelin biosynthesis and similarity (S), % of identical matches between reference PAO1 protein and its homolog, and length coverage (L), % of reference PAO1 protein sequence covered by its homolog. Gray box indicates there is no homolog under e-value 0.01 in the blastp search.

https://doi.org/10.1371/journal.ppat.1009681.s006

(TIF)

S1 Data. Information of 636 genomes.

https://doi.org/10.1371/journal.ppat.1009681.s007

(XLSX)

S2 Data. Pyseer results.

https://doi.org/10.1371/journal.ppat.1009681.s008

(XLSX)

S3 Data. Cluster and locus information of 636 genomes.

https://doi.org/10.1371/journal.ppat.1009681.s009

(XLSX)

S4 Data. Amino acid sequences of small insertion variants from isolates other than PAO1.

https://doi.org/10.1371/journal.ppat.1009681.s010

(DOCX)

S5 Data. Blastp results with 636 genomes and 238 YdjN homologs.

https://doi.org/10.1371/journal.ppat.1009681.s011

(XLSX)

S1 Table. Supplementary tables.

https://doi.org/10.1371/journal.ppat.1009681.s012

(DOCX)

Citation: Hwang W, Yong JH, Min KB, Lee K-M, Pascoe B, Sheppard SK, et al. (2021) Genome-wide association study of signature genetic alterations among pseudomonas aeruginosa cystic fibrosis isolates. PLoS Pathog 17(6): e1009681. https://doi.org/10.1371/journal.ppat.1009681

References

1. Stover CK, Pham XQ, Erwin AL, Mizoguchi SD, Warrener P, Hickey MJ, et al. Complete genome sequence of Pseudomonas aeruginosa PAO1, an opportunistic pathogen. Nature. 2000;406(6799):959–64. pmid:10984043

2. Bhagirath AY, Li Y, Somayajula D, Dadashi M, Badr S, Duan K. Cystic fibrosis lung environment and Pseudomonas aeruginosa infection. BMC Pulm Med. 2016;16(1):174-. pmid:27919253.

3. Silby MW, Winstanley C, Godfrey SAC, Levy SB, Jackson RW. Pseudomonas genomes: diverse and adaptable. FEMS Microbiology Reviews. 2011;35(4):652–80. pmid:21361996

4. Smith EE, Buckley DG, Wu Z, Saenphimmachak C, Hoffman LR, D’Argenio DA, et al. Genetic adaptation by Pseudomonas aeruginosa to the airways of cystic fibrosis patients. Proc Natl Acad Sci U S A. 2006;103(22):8487–92. Epub 2006/05/10. pmid:16687478.

5. Palmer KL, Mashburn LM, Singh PK, Whiteley M. Cystic fibrosis sputum supports growth and cues key aspects of Pseudomonas aeruginosa physiology. J Bacteriol. 2005;187(15):5267–77. pmid:16030221.

6. Fung C, Naughton S, Turnbull L, Tingpej P, Rose B, Arthur J, et al. Gene expression of Pseudomonas aeruginosa in a mucin-containing synthetic growth medium mimicking cystic fibrosis lung sputum. J Med Microbiol. 2010;59(Pt 9):1089–100. pmid:20522626.

7. Hoffman LR, Kulasekara HD, Emerson J, Houston LS, Burns JL, Ramsey BW, et al. Pseudomonas aeruginosa lasR mutants are associated with cystic fibrosis lung disease progression. J Cyst Fibros. 2009;8(1):66–70. Epub 2008/11/01. pmid:18974024.

8. Feltner JB, Wolter DJ, Pope CE, Groleau MC, Smalley NE, Greenberg EP, et al. LasR Variant Cystic Fibrosis Isolates Reveal an Adaptable Quorum-Sensing Hierarchy in Pseudomonas aeruginosa. mBio. 2016;7(5). Epub 2016/10/06. pmid:27703072.

9. Salunkhe P, Smart CH, Morgan JA, Panagea S, Walshaw MJ, Hart CA, et al. A cystic fibrosis epidemic strain of Pseudomonas aeruginosa displays enhanced virulence and antimicrobial resistance. J Bacteriol. 2005;187(14):4908–20. Epub 2005/07/05. pmid:15995206.

10. Mustafa MH, Khandekar S, Tunney MM, Elborn JS, Kahl BC, Denis O, et al. Acquired resistance to macrolides in Pseudomonas aeruginosa from cystic fibrosis patients. Eur Respir J. 2017;49(5). Epub 2017/05/21. pmid:28526799.

11. Maunders EA, Triniman RC, Western J, Rahman T, Welch M. Global reprogramming of virulence and antibiotic resistance in Pseudomonas aeruginosa by a single nucleotide polymorphism in elongation factor, fusA1. J Biol Chem. 2020. Epub 2020/09/19. pmid:32943550.

12. Barth AL, Pitt TL. Auxotrophic variants of Pseudomonas aeruginosa are selected from prototrophic wild-type strains in respiratory infections in patients with cystic fibrosis. J Clin Microbiol. 1995;33(1):37–40. pmid:7699062.

13. Barth A, Pitt TL. The high amino-acid content of sputum from cystic fibrosis patients promotes growth of auxotrophic Pseudomonas aeruginosa. J Med Microbiol. 1996;45:110–9. pmid:8683546

14. Marvig RL, Sommer LM, Molin S, Johansen HK. Convergent evolution and adaptation of Pseudomonas aeruginosa within patients with cystic fibrosis. Nature Genetics. 2015;47(1):57–64. pmid:25401299

15. Quinn RA, Lim YW, Mak TD, Whiteson K, Furlan M, Conrad D, et al. Metabolomics of pulmonary exacerbations reveals the personalized nature of cystic fibrosis disease. PeerJ. 2016;4:e2174–e. pmid:27602256.

16. Diaz Caballero J, Clark ST, Coburn B, Zhang Y, Wang PW, Donaldson SL, et al. Selective Sweeps and Parallel Pathoadaptation Drive Pseudomonas aeruginosa Evolution in the Cystic Fibrosis Lung. mBio. 2015;6(5):e00981. pmid:26330513.

17. La Rosa R, Johansen HK, Molin S. Convergent Metabolic Specialization through Distinct Evolutionary Paths in Pseudomonas aeruginosa. mBio. 2018;9(2):e00269–18. pmid:29636437.

18. Rosenberg NA, Huang L, Jewett EM, Szpiech ZA, Jankovic I, Boehnke M. Genome-wide association studies in diverse populations. Nat Rev Genet. 2010;11(5):356–66. Epub 2010/04/17. pmid:20395969.

19. Buniello A, MacArthur JAL, Cerezo M, Harris LW, Hayhurst J, Malangone C, et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 2019;47(D1):D1005–D12. pmid:30445434.

20. Hwang W, Yoon SS. Virulence Characteristics and an Action Mode of Antibiotic Resistance in Multidrug-Resistant Pseudomonas aeruginosa. Sci Rep. 2019;9(1):487. Epub 2019/01/27. pmid:30679735.

21. Chen PE, Shapiro BJ. The advent of genome-wide association studies for bacteria. Current Opinion in Microbiology. 2015;25:17–24. pmid:25835153

22. Lees JA, Vehkala M, Välimäki N, Harris SR, Chewapreecha C, Croucher NJ, et al. Sequence element enrichment analysis to determine the genetic basis of bacterial phenotypes. Nature Communications. 2016;7(1):12797. pmid:27633831

23. Sheppard SK, Didelot X, Meric G, Torralbo A, Jolley KA, Kelly DJ, et al. Genome-wide association study identifies vitamin B5 biosynthesis as a host specificity factor in Campylobacter. Proc Natl Acad Sci U S A. 2013;110(29):11923–7. Epub 2013/07/01. pmid:23818615.

24. Sheppard SK, Guttman DS, Fitzgerald JR. Population genomics of bacterial host adaptation. Nature Reviews Genetics. 2018;19(9):549–65. pmid:29973680

25. Earle SG, Wu C-H, Charlesworth J, Stoesser N, Gordon NC, Walker TM, et al. Identifying lineage effects when controlling for population structure improves power in bacterial association studies. Nat Microbiol. 2016;1:16041-. pmid:27572646.

26. Jaillard M, Lima L, Tournoud M, Mahe P, van Belkum A, Lacroix V, et al. A fast and agnostic method for bacterial genome-wide association studies: Bridging the gap between k-mers and genetic events. PLoS Genet. 2018;14(11):e1007758. Epub 2018/11/13. pmid:30419019 following competing interests: MJ, MT, PM and AvB are employees of bioMerieux, a company that develops and sells diagnostic tests in the field of infectious diseases. However, the study was designed and executed in an open manner and the presented method as well as all data generated have been deposited in the public domain, also resulting in the current publication.

27. Winsor GL, Griffiths EJ, Lo R, Dhillon BK, Shay JA, Brinkman FSL. Enhanced annotations and features for comparing thousands of Pseudomonas genomes in the Pseudomonas genome database. Nucleic Acids Res. 2016;44(D1):D646–D53. Epub 2015/11/17. pmid:26578582.

28. Ozer EA, Nnah E, Didelot X, Whitaker RJ, Hauser AR. The Population Structure of Pseudomonas aeruginosa Is Characterized by Genetic Isolation of exoU+ and exoS+ Lineages. Genome Biol Evol. 2019;11(1):1780–96. pmid:31173069.

29. Lees JA, Galardini M, Bentley SD, Weiser JN, Corander J. pyseer: a comprehensive tool for microbial pangenome-wide association studies. Bioinformatics. 2018;34(24):4310–2. pmid:30535304

30. Hadfield J, Croucher NJ, Goater RJ, Abudahab K, Aanensen DM, Harris SR. Phandango: an interactive viewer for bacterial population genomics. Bioinformatics. 2018;34(2):292–3. pmid:29028899

31. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30(14):2068–9. pmid:24642063

32. Anthony M, Rose B, Pegler MB, Elkins M, Service H, Thamotharampillai K, et al. Genetic analysis of Pseudomonas aeruginosa isolates from the sputa of Australian adult cystic fibrosis patients. J Clin Microbiol. 2002;40(8):2772–8. pmid:12149328.

33. Boucher JC, Yu H, Mudd M, Deretic V. Mucoid Pseudomonas aeruginosa in cystic fibrosis: Characterization of muc mutations in clinical isolates and analysis of clearance in a mouse model of respiratory infection. Infection and immunity. 1997;65:3838–46. pmid:9284161

34. Kanehisa M, Sato Y, Morishima K. BlastKOALA and GhostKOALA: KEGG Tools for Functional Characterization of Genome and Metagenome Sequences. Journal of Molecular Biology. 2016;428(4):726–31. pmid:26585406

35. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. pmid:10592173.

36. Yoon SS, Coakley R, Lau GW, Lymar SV, Gaston B, Karabulut AC, et al. Anaerobic killing of mucoid Pseudomonas aeruginosa by acidified nitrite derivatives under cystic fibrosis airway conditions. J Clin Invest. 2006;116(2):436–46. Epub 2006/01/26. pmid:16440061.

37. Min KB, Lee K-M, Oh YT, Yoon SS. Nonmucoid conversion of mucoid Pseudomonas aeruginosa induced by sulfate-stimulated growth. FEMS Microbiology Letters. 2014;360(2):157–66. pmid:25227776

38. Morita Y, Tomida J, Kawamura Y. MexXY multidrug efflux system of Pseudomonas aeruginosa. Front Microbiol. 2012;3:408-. pmid:23233851.

39. Rojo F. Carbon catabolite repression in Pseudomonas: optimizing metabolic versatility and interactions with the environment. FEMS Microbiology Reviews. 2010;34(5):658–84. pmid:20412307

40. Barraud N, Buson A, Jarolimek W, Rice SA. Mannitol enhances antibiotic sensitivity of persister bacteria in Pseudomonas aeruginosa biofilms. PLoS One. 2013;8(12):e84220–e. pmid:24349568.

41. Gough J, Karplus K, Hughey R, Chothia C. Assignment of homology to genome sequences using a library of hidden Markov models that represent all proteins of known structure11Edited by G. Von Heijne. Journal of Molecular Biology. 2001;313(4):903–19. pmid:11697912

42. Dolan SK, Pereira G, Silva-Rocha R, Welch M. Transcriptional regulation of central carbon metabolism in Pseudomonas aeruginosa. Microb Biotechnol. 2020;13(1):285–9. Epub 2019/06/11. pmid:31187593.

43. Campilongo R, Fung RKY, Little RH, Grenga L, Trampari E, Pepe S, et al. One ligand, two regulators and three binding sites: How KDPG controls primary carbon metabolism in Pseudomonas. PLOS Genetics. 2017;13(6):e1006839. pmid:28658302

44. Ciofu O, Tolker-Nielsen T. Tolerance and Resistance of Pseudomonas aeruginosa Biofilms to Antimicrobial Agents-How P. aeruginosa Can Escape Antibiotics. Front Microbiol. 2019;10:913-. pmid:31130925.

45. Liu GY, Nizet V. Color me bad: microbial pigments as virulence factors. Trends Microbiol. 2009;17(9):406–13. Epub 2009/08/31. pmid:19726196.

46. Cox CD. Role of pyocyanin in the acquisition of iron from transferrin. Infection and immunity. 1986;52(1):263–70. pmid:2937736.

47. Ghio AJ, Roggli VL, Soukup JM, Richards JH, Randell SH, Muhlebach MS. Iron accumulates in the lavage and explanted lungs of cystic fibrosis patients. Journal of Cystic Fibrosis. 2013;12(4):390–8. pmid:23176785

48. Meyer KC, Sharma A, Brown R, Weatherly M, Moya FR, Lewandoski J, et al. Function and Composition of Pulmonary Surfactant and Surfactant-Derived Fatty Acid Profiles Are Altered in Young Adults With Cystic Fibrosis. Chest. 2000;118(1):164–74. pmid:10893374

49. Renilla S, Bernal V, Fuhrer T, Castaño-Cerezo S, Pastor JM, Iborra JL, et al. Acetate scavenging activity in Escherichia coli: interplay of acetyl–CoA synthetase and the PEP–glyoxylate cycle in chemostat cultures. Applied Microbiology and Biotechnology. 2012;93(5):2109–24. pmid:21881893

50. Son MS, Matthews WJ Jr, Kang Y, Nguyen DT, Hoang TT. In vivo evidence of Pseudomonas aeruginosa nutrient acquisition and pathogenesis in the lungs of cystic fibrosis patients. Infection and immunity. 2007;75(11):5313–24. Epub 2007/08/27. pmid:17724070.

51. Flynn JM, Phan C, Hunter RC. Genome-Wide Survey of Pseudomonas aeruginosa PA14 Reveals a Role for the Glyoxylate Pathway and Extracellular Proteases in the Utilization of Mucin. Infection and immunity. 2017;85(8):e00182–17. pmid:28507068.

52. Rossi E, Falcone M, Molin S, Johansen HK. High-resolution in situ transcriptomics of Pseudomonas aeruginosa unveils genotype independent patho-phenotypes in cystic fibrosis lungs. Nature Communications. 2018;9(1):3459. pmid:30150613

53. Lorenz MC, Fink GR. Life and death in a macrophage: role of the glyoxylate cycle in virulence. Eukaryot Cell. 2002;1(5):657–62. pmid:12455685.

54. Lemire J, Alhasawi A, Appanna VP, Tharmalingam S, Appanna VD. Metabolic defence against oxidative stress: the road less travelled so far. Journal of Applied Microbiology. 2017;123(4):798–809. pmid:28609580

55. Donlan RM. Biofilms: microbial life on surfaces. Emerg Infect Dis. 2002;8(9):881–90. pmid:12194761.

56. Käll L, Krogh A, Sonnhammer ELL. Advantages of combined transmembrane topology and signal peptide prediction—the Phobius web server. Nucleic Acids Res. 2007;35(Web Server issue):W429–W32. Epub 2007/05/05. pmid:17483518.

57. Ohtsu I, Kawano Y, Suzuki M, Morigasaki S, Saiki K, Yamazaki S, et al. Uptake of L-cystine via an ABC transporter contributes defense of oxidative stress in the L-cystine export-dependent manner in Escherichia coli. PLoS One. 2015;10(3):e0120619–e. pmid:25837721.

58. The UniProt C. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47(D1):D506–D15. pmid:30395287

59. Palmer KL, Aye LM, Whiteley M. Nutritional cues control Pseudomonas aeruginosa multicellular behavior in cystic fibrosis sputum. J Bacteriol. 2007;189(22):8079–87. Epub 2007/09/14. pmid:17873029.

60. Korshunov S, Imlay KRC, Imlay JA. Cystine import is a valuable but risky process whose hazards Escherichia coli minimizes by inducing a cysteine exporter. Molecular Microbiology. 2020;113(1):22–39. pmid:31612555

61. Cornelis P, Dingemans J. Pseudomonas aeruginosa adapts its iron uptake strategies in function of the type of infections. Front Cell Infect Microbiol. 2013;3:75-. pmid:24294593.

62. Ronnebaum TA, Lamb AL. Nonribosomal peptides for iron acquisition: pyochelin biosynthesis as a case study. Current Opinion in Structural Biology. 2018;53:1–11. pmid:29455106

63. Lee K, Lee K-M, Go J, Ryu J-C, Ryu J-H, Yoon SS. The ferrichrome receptor A as a new target for Pseudomonas aeruginosa virulence attenuation. FEMS Microbiology Letters. 2016;363(11). pmid:27190289

64. Ochsner U, Johnson Z, Vasil M. Genetics and regulation of two distinct haem-uptake systems, phu and has, in Pseudomonas aeruginosa. Microbiology (Reading, England). 2000;146 (Pt 1):185–98. pmid:10658665

65. Marvig RL, Damkiær S, Khademi SMH, Markussen TM, Molin S, Jelsbak L. Within-host evolution of Pseudomonas aeruginosa reveals adaptation toward iron acquisition from hemoglobin. mBio. 2014;5(3):e00966. pmid:24803516.

66. Kobras CM, Fenton AK, Sheppard SK. Next-generation microbiology: from comparative genomics to gene function. Genome Biol. 2021;22(1):123-. pmid:33926534.

67. Simonsen M, Mailund T, Pedersen CNS, editors. Rapid Neighbour-Joining. Algorithms in Bioinformatics; 2008 2008//; Berlin, Heidelberg: Springer Berlin Heidelberg.

68. Seemann T. Snippy. https://github.com/tseemann/snippy.

69. Argimón S, Abudahab K, Goater RJE, Fedosejev A, Bhai J, Glasner C, et al. Microreact: visualizing and sharing data for genomic epidemiology and phylogeography. Microb Genom. 2016;2(11):e000093–e. pmid:28348833.

70. Menardo F, Loiseau C, Brites D, Coscolla M, Gygli SM, Rutaihwa LK, et al. Treemmer: a tool to reduce large phylogenetic datasets with minimal loss of diversity. BMC Bioinformatics. 2018;19(1):164. pmid:29716518

71. Välimäki N. Fsm-lite. https://github.com/nvalimak/fsm-lite.

72. Bayliss S. Snp-sim. https://github.com/SionBayliss/bio-scripts/tree/master/scripts.

73. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52. pmid:21572440.

74. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC bioinformatics. 2009;10:421-. pmid:20003500.

75. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80. Epub 2013/01/16. pmid:23329690.

76. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994;22(22):4673–80. pmid:7984417.

77. Waterhouse AM, Procter JB, Martin DMA, Clamp M, Barton GJ. Jalview Version 2—a multiple sequence alignment editor and analysis workbench. Bioinformatics (Oxford, England). 2009;25(9):1189–91. Epub 2009/01/16. pmid:19151095.

78. Spilker T, LiPuma JJ. Draft Genome Sequences of 63 Pseudomonas aeruginosa Isolates Recovered from Cystic Fibrosis Sputum. Genome Announc. 2016;4(2):e00231–16. pmid:27103710.

79. Choi K-H, Gaynor JB, White KG, Lopez C, Bosio CM, Karkhoff-Schweizer RR, et al. A Tn7-based broad-range bacterial cloning and expression system. Nature Methods. 2005;2(6):443–8. pmid:15908923

80. Sezgi S, Seyhan U, Gulgun B-T, Aysegul Y. Quorum sensing and virulence of Pseudomonas aeruginosa during urinary tract infections. The Journal of Infection in Developing Countries. 2012;6(06). pmid:22706193

81. Carlsson M, Shukla S, Petersson AC, Segelmark M, Hellmark T. Pseudomonas aeruginosa in cystic fibrosis: Pyocyanin negative strains are associated with BPI-ANCA and progressive lung disease. Journal of Cystic Fibrosis. 2011;10(4):265–71. pmid:21463973

82. Kim E-J, Wang W, Deckwer W-D, Zeng A-P. Expression of the quorum-sensing regulatory protein LasR is strongly affected by iron and oxygen concentrations in cultures of Pseudomonas aeruginosa irrespective of cell density. Microbiology (Reading, England). 2005;151:1127–38. pmid:15817780

83. O’Toole GA. Microtiter dish biofilm formation assay. J Vis Exp. 2011;(47):2437. pmid:21307833.

84. Münch R, Hiller K, Barg H, Heldt D, Linz S, Wingender E, et al. PRODORIC: prokaryotic database of gene regulation. Nucleic Acids Res. 2003;31(1):266–9. pmid:12519998

85. Crooks GE, Hon G, Chandonia J-M, Brenner SE. WebLogo: a sequence logo generator. Genome Res. 2004;14(6):1188–90. pmid:15173120.