Parallel Evolution in Pseudomonas aeruginosa over 39,000 Generations In Vivo

Abstract

Introduction Pathogenic microbes and viruses, including those that cause AIDS, tuberculosis, and malaria, can remain closely associated with their hosts for decades. While much has been learned about microbial evolution from long-term in vitro experiments (reviewed in reference 1), less is known about pathogen evolution during chronic infection. Pseudomonas aeruginosa, a ubiquitous Gram-negative bacterium, is a common cause of chronic illness in individuals with the heritable disease cystic fibrosis (CF), where it is the leading cause of morbidity and mortality (2). Upon infection, individuals with CF are unable to clear P. aeruginosa from the lungs, where it grows to high cell densities despite competing microorganisms, a strong host immune response, and antibiotic treatment (2). Several important features make P. aeruginosa CF lung infections suitable for studying microbial evolution in vivo. First, a single clone can colonize the lungs of a patient with CF for over a decade (3), enabling the study of chronological clonal isolates. In addition, much is known about the CF lung environment, including carbon source availability, oxygen content, and in vivo growth rates (4–10). Previous studies have proposed that P. aeruginosa undergoes adaptive evolution in the CF lung primarily based on two observations. (i) There is a high ratio of nonsynonymous to synonymous mutations that occur over time (11). (ii) Isolates from chronic infections evolve similar phenotypes, including the following: formation of small colonies on agar (small colony variants or dwarf strains) (12, 13), overproduction of alginate exopolysaccharide (mucoid strains) (12, 14), quorum-sensing deficiency (11, 15), motility reduction (16, 17), altered lipopolysaccharide (18), reduced virulence factor expression (19, 20), and hypermutation (21). However, it is unclear how global changes in gene expression mediate these “CF-evolved” phenotypes. Transcriptional profiling has been used to identify potential adaptive traits in two similarly grown Escherichia coli lineages after 20,000 generations (22). Based on these in vitro experiments, we hypothesized that transcriptional profiling would identify potential adaptive expression traits in P. aeruginosa lineages that evolved in vivo. As detailed analysis of gene expression profiling has not been performed on sequential P. aeruginosa strains obtained from multiple individuals with CF, we conducted transcriptomic analyses on 17 chronological clonal isolates collected over 3 months to 8 years from 3 CF patients. This reflects approximately 1,200 to 39,000 generations based on a previously estimated P. aeruginosa doubling time of 100 min in the CF lung (10). Our studies reveal that parallel changes in gene expression occur during in vivo evolution and that P. aeruginosa uses multiple pathways to establish chronic infection. RESULTS In order to identify gene expression changes that occur during chronic infection, transcriptional profiling was carried out on clonal P. aeruginosa isolates collected from 3 patients with CF (Fig. 1). We anticipated that these analyses would reveal how different molecular processes are coordinated to establish long-term infection. Clonal isolates were collected between 3 months and 8 years after colonization. Individuals A, B, and C harbored unique P. aeruginosa clones (data not shown). A strain replacement occurred in individual C between 1983 and 1987, so these isolates are divided into clonal groups Ca and Cb. FIG 1 Pseudomonas aeruginosa isolates from three individuals with cystic fibrosis. Note that clonal groups Ca and Cb are isolated from the same individual, though they are different clonal groups (confirmed by RAPD assay and pulsed-field gel electrophoresis [data not shown]). More information is available in Table S1 in the supplemental material. SCV, small colony variant. Microarray analyses. The Affymetrix P. aeruginosa GeneChip was designed from the P. aeruginosa PAO1 genome; therefore, it may be argued that the changes we observed reflect genetic variations between isolates rather than gene expression changes. However, a number of factors do not support this claim. First, the Affymetrix P. aeruginosa GeneChip is designed to capture genes using 13 unique 25-mer probes. Therefore, small variations (such as single-nucleotide polymorphisms) would not be significant after probe summarization. Second, we considered only genes that have orthologs in four other P. aeruginosa genomes (see Materials and Methods), so genes that are variable between P. aeruginosa strains were not included in our analyses. Third, the same microarray platform has been used successfully to assess the genomic content of various P. aeruginosa isolates, including environmental and CF isolates (23). Based on genomic DNA hybridization to the Affymetrix P. aeruginosa PAO1 GeneChip, it was reported that between 96.1% and 97.7% of the P. aeruginosa genome is conserved, including almost all known virulence factors (23). These data suggest that the Affymetrix P. aeruginosa PAO1 GeneChip is able to capture the conserved region of each gene in most P. aeruginosa strains. Finally, we analyzed our data using the Affymetrix MAS5 algorithm, which calls the presence or absence of target cDNAs based on the distribution of perfect match and mismatch probe signals. If our observations were derived from mutations that accumulated over time, we would have observed a high incidence of absence calls in the late strains compared to the early strains. However, as we report in Table S6 in the supplemental material, we did not see this trend. These results clearly show that our observations reflect changes in gene expression rather than genetic variations. Clonal isolates cluster by patient. To view the global features of our transcriptomic data, we first conducted hierarchical clustering of microarray signals obtained from each isolate. We used two similarity measures, Spearman correlation coefficient and Euclidean distance, to overcome the bias that can occur with one similarity measurement. Overall, transcriptomic clustering (Fig. 2) clearly showed that isolates clustered by clonal group (P = 0.001 by analysis of similarity [ANOSIM]) rather than by morphological phenotype (P = 0.263 by ANOSIM) or time in the CF lung (P = 0.09 by ANOSIM). Similar results were obtained with both clustering methods. These data indicate that despite growth in the lungs for thousands of generations, P. aeruginosa isolates collected at later time points, at least in regard to gene expression, resemble the ancestor more than isolates from other individuals. FIG 2 Hierarchical clustering of microarray data. Normalized raw microarray signals of 5,391 genes (mean values for two biological replicates) were used for clustering. We confirmed that clusters of each clonal group are not dependent on the similarity measurement method by using both a Spearman correlation coefficient (A) and Euclidean distance (B). Multiple solutions for establishing chronic infection exist. Transcriptional profiling showed that clonal groups are highly similar, but there are some differences that account for each group clustering by patient. In order to identify these differences, we compared clonal groups (including ancestor strains) to each other. Clonal group A differentially expressed 65 genes, clonal group B differentially expressed 58 genes, clonal group Ca differentially expressed 119 genes, and clonal group Cb differentially expressed 99 genes (see Table S4 in the supplemental material). A majority of these genes were involved in virulence, quorum sensing, alginate production, branched-chain amino acid metabolism, and motility (Fig. 3; see Table S4 in the supplemental material). These data suggest that in regard to gene expression, there are multiple solutions for establishing chronic infection in the lungs of patients with CF. FIG 3 Microarray heat maps of genes differentially expressed between clonal groups. Differentially expressed genes between clonal group A (A), clonal group B (B), clonal group Ca (C), and clonal group Cb (D) are shown. Each heat map shows a subset of differentially expressed genes between clonal groups with signals of >1,000 in at least one strain. Comprehensive lists for these analyses are available in Table S4 in the supplemental material. Red indicates high levels of mRNA, and blue indicates low levels of mRNA. Ancestor strains are early colonizers of the CF lung. To our knowledge, the ancestors were the first chronically established strains in each individual. However, these strains may have undergone adaptation before they were collected. In order to show that ancestral strains had not already undergone significant adaptation within the CF lung, we compared ancestor strains A1, B1, and Ca1 to the laboratory reference strain, PA14. PA14 was isolated from an acute burn wound; it has undergone minimal passage in the laboratory and maintains many of its original traits, making it a valuable non-CF reference strain. Ancestor strains A1, B1, and Ca1 differentially expressed 228, 181, and 265 genes, respectively, compared to strain PA14 (see Table S5 in the supplemental material). In contrast, late isolates showed more than 600 gene expression changes compared to strain PA14. These results suggest that the ancestral strains show moderate gene expression changes compared to strain PA14, but since later isolates showed approximately 3-fold-more differences, the strains are still adapting to the CF lung. Strain Cb1 exhibited an adaptive phenotype, mucoidy, and is likely not the original ancestor of subsequent isolates, although this strain does display adaptation (see below). Clonal groups exhibit parallel changes in gene expression. We next examined gene expression changes within each lineage by comparing the transcriptomes of all isolates in each clonal group to their ancestor (Fig. 4; see Table S6 in the supplemental material). We focused on differentially expressed genes that were maintained throughout infection compared to each ancestor. Isolates in clonal group A differentially expressed 76 genes; these changes were maintained over 3 years. Four hundred forty-one changes were maintained over 8 years in isolates from clonal group B. Clonal group Ca differentially regulated 37 genes over 3 months, while clonal group Cb differentially regulated 281 genes over 7 years. These results indicate that 1 to 8% of P. aeruginosa genes displayed differences in expression after chronic infection within the CF lung. FIG 4 Microarray heat map of differentially expressed genes within clonal groups. All strains in each clonal group were compared to their ancestor (initial strain) and evaluated for significant changes in gene expression (>2-fold change and FDR of <0.05). The numbers used to generate this heat map are the mean fold changes (log2) in mRNA levels compared to the mRNA level of the ancestor strain of each clonal group. Red indicates an increase in the mRNA level, and blue indicates a decrease in the mRNA level. All genes detected in this analysis are reported separately in Table S6 in the supplemental material. Although the environment of the CF lung likely differs between individuals, we hypothesized that due to strong selective pressures in all patients (immune response, treatment, and consistent nutritional environment), P. aeruginosa isolates from different patients would display common changes in gene expression. To test this hypothesis, we compared the above within-group analyses and identified a set of genes that were commonly differentially expressed within all 3 lineages (Table 1). Clonal group Ca evolved over a short time period, 3 months; therefore, it was excluded from the analysis. Twenty-four genes commonly changed over time in clonal groups A, B, and Cb. Importantly, these genes changed in the same direction across all 3 lineages. Nine genes were downregulated, and 15 were upregulated. Four of the downregulated genes encoded proteins of unknown function, while the remaining were involved in type 4 fimbrial biogenesis. Ten genes encoding proteins of unknown function were upregulated. The other upregulated genes encode two outer membrane proteins (PA1048, OsmE), PA4880 (probable bacterioferritin), phaF (polyhydroxyalkanoate synthesis protein PhaF), and PA1562 (aconitase). Additionally, 86 genes were commonly regulated in at least 2 lineages (Fig. 5; see Table S6 in the supplemental material). Many of these genes were involved in flagellar biosynthesis, fimbrial (pilus) biosynthesis, and polyamine transport. The appearance of similar alterations in gene expression patterns, particularly changes in the same direction, suggests that P. aeruginosa undergoes parallel evolution in the CF lung. FIG 5 Common gene expression changes among three clonal groups. Venn diagram showing common changes in gene expression between three clonal groups (groups A, B, and Cb). The 24 genes commonly regulated in the three clonal groups are listed in Table 1. Genes that are positively regulated over time (+) and genes that are negatively regulated over time (-) are indicated (for instance, in the red circle at the top of the figure, of the 16 genes, 4 genes were positively regulated over time and 12 genes were negatively regulated over time.) The P value for the two clonal groups was calculated using a hypergeometric test and indicates that the probability of identifying the overlapping genes by random chance is very low. The full list of genes is available in Table S6 in the supplemental material. TABLE 1 Genes exhibiting parallel expression changes in three clonal groups a Gene b Annotation b Change in expression c PA0045 Hypothetical protein Decrease PA0046 Hypothetical protein Decrease PA0047 Hypothetical protein Decrease PA0411 (pilJ) Type 4 fimbrial biogenesis protein Decrease PA5041 (pilP) Type 4 fimbrial biogenesis protein Decrease PA5042 (pilO) Type 4 fimbrial biogenesis protein Decrease PA5043 (pilN) Type 4 fimbrial biogenesis protein Decrease PA5044 (pilM) Type 4 fimbrial biogenesis protein Decrease PA5139 Hypothetical protein Decrease PA1048 Probable outer membrane protein precursor Increase PA1106 Hypothetical protein Increase PA1323 Hypothetical protein Increase PA1324 Hypothetical protein Increase PA1471 Hypothetical protein Increase PA1562 (acnA) Aconitase Increase PA1592 Hypothetical protein Increase PA2485 Hypothetical protein Increase PA2779 Hypothetical protein Increase PA3040 Conserved hypothetical protein Increase PA3691 Hypothetical protein Increase PA4876 (osmE) Osmotically inducible lipoprotein Increase PA4880 Probable bacterioferritin Increase PA5060 (phaF) Polyhydroxyalkanoate synthesis protein Increase PA5178 Conserved hypothetical protein Increase a  Only genes that have a >2-fold change in mRNA levels compared to the mRNA level of each ancestor strain (FDR of <0.05), changed in clonal groups A, B, and Cb, and maintained changes over time are listed. The raw data are available in Table S6 in the supplemental material. b  Annotation data were downloaded from the Pseudomonas Genome Database on 23 November 2009. c  mRNA levels of the late strains were compared to those of the ancestor strain in each clonal group. “Increase” indicates that mRNA levels were higher in the late strains than in the ancestor strain, and “decrease” indicates that mRNA levels were lower in the late strains than in the ancestor strain. DISCUSSION Our goal was to identify adaptive expression traits of P. aeruginosa during chronic CF lung infection. Our study is limited by the absence of genetic data because we cannot correlate specific mutations with changes in gene expression, and pleiotropic effects cannot be ruled out. However, since genetic mutations can affect gene expression, our study is important for understanding the global evolution of P. aeruginosa during chronic infection. Our transcriptional analysis is the most comprehensive study of this type thus far. While D’Argenio et al. (24) and Hoboth et al. (25) compared transcriptional profiles of early and late isolates, these studies were limited by the number of samples or patients. In contrast, we characterized chronological, clonal isolates collected from multiple CF patients. We identified 24 genes that showed similar changes in gene expression across 3 separate P. aeruginosa lineages. This is striking because our statistical analyses show that overlap between groups is highly significant (P value ≤10−10) (Fig. 5), so the probability of identifying 24 genes by random chance is very low. Since patients with CF can carry diverse populations of P. aeruginosa at any given time, one obvious limitation to this study is sampling. However, all 3 late isolates recovered on the same date from patient B (B3.1, B3.2, and B3.3) showed 23 of 24 common gene expression changes, which strongly supports parallel evolution. Parallel evolution occurs when two closely related organisms independently develop the same adaptive traits due to the nature of their environments (reviewed in reference 26). Parallelism is indicative of adaptive evolution, which has been proposed to occur in the CF lung (11). Our data are similar to those of Cooper et al. (22), who demonstrated that 2 independently evolving E. coli populations growing under similar laboratory conditions displayed the same 59 changes in gene expression after 20,000 generations. Later, Barrick et al. (27) showed that most of the 45 mutations occurring in this long-term evolution experiment were beneficial. On the basis of these data, we hypothesize that some of the 24 commonly regulated genes encode adaptive traits (Table 1). There are several lines of evidence to support this claim. First, strong parallelism is a good indicator of adaptive evolution. Second, 5 fimbrial (pilus) biosynthetic genes (pilJ, pilP, pilO, pilN, and pilM) were downregulated over time in all 3 lineages (Table 1). It is well documented that the loss of pili protects P. aeruginosa from phagocytosis in vivo by neutrophils, thus allowing P. aeruginosa to escape a primary immune component in the lungs of patients with CF (16, 28–30). Third, 4 genes important for P. aeruginosa biofilm formation in vitro (PA5139, PA1592, PA2779, and PA4876) and 3 genes upregulated during biofilm growth in vitro (PA1471, PA4876, and PA1324) were differentially regulated over time in all 3 lineages (31–33). Resistance to antibiotics and to host clearance is partly attributed to P. aeruginosa biofilm growth in the CF lung (34–36), suggesting that adaptation to the biofilm lifestyle via differential expression of these genes is important for maintaining long-term infections. Finally, of the 24 genes, 3 (PA4876, PA1323, and PA1324) were upregulated in an important P. aeruginosa CF isolate, the Liverpool epidemic strain (LES), compared to laboratory strain PAO1 (37). LES is a particularly virulent P. aeruginosa CF strain; thus, expression of these traits may represent an important adaptation to the lung environment. Our results indicate that parallel evolution occurs in the lungs of individuals with CF; however, our data also show that specific, within-lineage changes also occur. The contribution of parallel changes versus within-lineage changes is unclear. The number of traits identified as evolving in parallel (24) is smaller than the number of changes in gene expression within lineages (up to 441) and may represent only a core set of changes (Fig. 5). For example, distinct virulence traits appear to be utilized by different clonal groups to establish and maintain infection, supporting previous hypotheses that P. aeruginosa employs unique evolutionary pathways to establish chronic infection (38) and that genes required for pathogenicity in one strain may not be required or predictive for others (38–40). Thus, we propose that natural selection may not be acting on individual virulence traits, but instead on some or all of the 24 genes identified in our study. In microbes, parallel gene expression changes have been observed only in long-term in vitro evolution experiments. Our study shows parallel changes in gene expression of strains grown in vivo for up to 39,000 generations. It should be noted that we measured gene expression in vitro. While it is well documented that the medium used to grow bacteria in this study (synthetic cystic fibrosis sputum medium [SCFM]) closely mimics the nutritional environment of the CF lung (7), the number of traits undergoing parallel evolution may be underestimated, or some may be expressed only under our culture conditions. It would be ideal to sample P. aeruginosa RNA directly from sputum samples from patients with CF, but the extraction process is challenging and would likely reflect species heterogeneity. Finally, 19 genes (of 24) identified in this study represent traits that have not been previously identified as adaptive in the CF lung and are therefore temporal markers for P. aeruginosa chronic colonization. Many of these genes encode proteins of unknown function, suggesting that future studies aimed at understanding the roles of these genes could provide unique insight into selective pressures in the CF lung. MATERIALS AND METHODS Isolate collection and genetic analyses. Clonal isolates of P. aeruginosa were collected and clonally purified from three patients with CF at the British Columbia’s Children’s Hospital, Shaughnessy Hospital, or St. Paul’s Hospital (all three hospitals are located in Vancouver, British Columbia, Canada) as previously described (16). Briefly, P. aeruginosa was isolated from sputum or throat samples and plated on Columbia agar supplemented with 5% sheep blood, MacConkey agar, or chocolate agar. On the basis of morphology, isolates were characterized as classic (C), dwarf (D), or mucoid (M) (12). Isolates A1, B1, and Ca2 were collected from respiratory sites, and the remaining isolates were collected from the throat or sputum (Fig. 1; see Table S1 in the supplemental material). The isolates were subcultured on Columbia blood agar, resuspended in 2 ml Mueller-Hinton broth with 8% dimethyl sulfoxide (DMSO), and stored at −75°C. To confirm that the isolates were P. aeruginosa, they were plated on a selective medium, FC agar (41). Genomic DNA extraction and genetic typing were carried out as previously described (3). P. aeruginosa isolates were typed by randomly amplified polymorphic DNA (RAPD) typing using primer 272 (3). All isolates from each individual are clonal except for patient C; a strain replacement occurred in patient C between 1983 and 1987 (Fig. 1; see Table S1 in the supplemental material). For clarity, clonal isolates from individuals A, B, and C will be referred to as clonal groups A, B, Ca, and Cb. Isolates collected from patient C are split into 2 clonal groups, Ca and Cb, due to the aforementioned strain replacement. Chronological isolates from each clonal set are designated by the patient letter (A, B, Ca, and Cb) along with the temporal order of isolation. For example, A1 was the first isolate collected from individual A, A2 was the second isolate, and A3.1 and A3.2 were isolated on the same date but displayed different morphologies (Fig. 1; see Table S1 in the supplemental material). P. aeruginosa annotation. All analyses were based on the P. aeruginosa PAO1 annotation (downloaded from the Pseudomonas Genome Database on 23 November 2009). To exclude PAO1 strain-specific genes, we defined a conserved gene set using 4 P. aeruginosa genomes maintained at PseudoCAP (PAO1, PA14, PA7, and LESB58) based on the reciprocal best-hit BLAST method. We identified 5,465 PAO1 genes (out of 5,569 PAO1 genes; 4,699 genes were conserved among all 4 strains) that have at least one ortholog in the other 3 strains; these genes were considered for further analysis (see Table S2 in the supplemental material). To assign microarray probe sets to these genes, we downloaded probe sequences from the Affymetrix website (Affymetrix, Santa Clara, CA) and mapped them to the P. aeruginosa PAO1 genome (GenBank accession no. NC002516.2) by using Exonerate (version 2.20) (42). After discarding probes that were not uniquely mapped on the genome, we mapped them again to P. aeruginosa PAO1 cDNA sequences. If less than 12 probes or more than 14 probes in a probe set were mapped to a gene, it was discarded, since each gene was represented by 13 unique probes. Thus, a total of 5,391 coding genes were considered in microarray data analysis (see Table S3 in the supplemental material). Expression profiling with Affymetrix microarrays. Microarray analyses were performed in duplicate on clonal isolates from 3 different individuals with CF. These analyses included the following strains: A1, A2, A3.1, A3.2, A4 (A strains from individual A), B1, B2.1, B2.2, B2.3, B3.1, B3.2, B3.3 (B strains from individual B), Ca1, Ca2, Cb1, Cb2, Cb3 (Ca and Cb strains from individual C), UCBPP-PA14 (PA14), and PAO1. Clinical isolates and the reference strains PA14 and PAO1 (6) were routinely grown on Difco blood agar base (BD Sciences) supplemented with 5% sheep blood (Remel) or brain heart infusion agar (Fisher). All Affymetrix GeneChip experiments were performed in synthetic cystic fibrosis sputum medium (SCFM), which mimics the nutritional environment of the CF lung (7). Bacterial growth in liquid medium (25 ml in a 250-ml flask) was monitored by measuring the optical density at 600 nm (OD600) during growth at 37°C with shaking at 250 rpm. Global gene expression profiling was carried out as previously described (6, 43) with minor modifications. P. aeruginosa isolates were grown in SCFM, and the cells were harvested during exponential phase (OD600 of 0.4 to 0.5). Cultures were mixed 1:1 with RNAlater (Ambion), an RNA-stabilizing agent. RNA was isolated using the RNeasy minikit (Qiagen), and cDNA was prepared for Affymetrix GeneChip microarray analysis as previously described (6, 43). PCR amplification of the P. aeruginosa rplU gene (44, 45) was used to detect DNA contamination using the primers rplU-For (named For for forward primer) (5′-CGCAGTGATTGTTACCGGTG-3′) and rplU-Rev (named Rev for reverse primer) (5′-AGGCCTGAATGCCGGTGATC-3′). To assess RNA integrity, the samples were subjected to agarose gel electrophoresis. Affymetrix GeneChips were washed, stained, and scanned using an Affymetrix fluidics station at the University of Iowa DNA core facility. Microarray analyses. We preprocessed microarray .CEL files by the RMA (robust multiarray analysis) method by using the affy package (version 1.18.2) (46) in R (version 2.8.1) with default options (correction with perfect match probes only, quantile normalization, and expression measure by median polish). Hierarchical clustering analysis was performed with these signals after summarizing two signals from biological replicates by their mean and tested by both Euclidean distance and Spearman correlation coefficient as similarity measures. Gene expression differences were evaluated by a linear model, implemented on the limma package in R (version 2.16.5) (47). If the maximum signal of the probe set among all isolates was less than 100, it was not considered in the analysis. We considered genes to be significantly differentially expressed if they changed >2-fold and the false discovery rate (FDR) was less than 0.05. Differentially expressed genes between clonal groups (see Table S4 in the supplemental material) were identified by the following 2 steps. (i) We compared the expression levels of all isolates in the same clonal group to those in another clonal group. (ii) We identified differentially expressed genes between two groups based on an FDR of <0.05 and a >2.0-fold change. Differentially expressed genes over time within clonal groups (see Table S6 in the supplemental material) were identified by the following five steps or criteria. (i) We compared the gene expression levels of all isolates from each clonal group to the ancestor. (ii) We identified differentially expressed genes based on an FDR of <0.05 and a >2.0-fold change. (iii) We discarded genes that changed in the direction opposite the direction of the ancestor (i.e., upregulated in the intermediate isolate and downregulated in the late isolate). (iv) If the level of gene expression was significantly different in any isolate within the same group, it was selected. (v) If a gene was expressed differently in both the intermediate and late isolates, it was selected. ANOSIM test was conducted by using the vegan (http://vegan.r-forge.r-project.org/) package in R (48). Microarray data accession number. All microarray data are available at the NCBI GEO database (accession no. GSE21966). Supplemental material is also available on the World Wide Web at http://www.marcottelab.org/index.php/PSEAE_CF.2010.