Comparative Population Genomics Analysis of the Mammalian Fungal Pathogen Pneumocystis

Abstract

Pneumocystis species are opportunistic mammalian pathogens that cause severe pneumonia in immunocompromised individuals. These fungi are highly host specific and uncultivable in vitro Human Pneumocystis infections present major challenges because of a limited therapeutic arsenal and the rise of drug resistance. To investigate the diversity and demographic history of natural populations of Pneumocystis infecting humans, rats, and mice, we performed whole-genome and large-scale multilocus sequencing of infected tissues collected in various geographic locations. Here, we detected reduced levels of recombination and variations in historical demography, which shape the global population structures. We report estimates of evolutionary rates, levels of genetic diversity, and population sizes. Molecular clock estimates indicate that Pneumocystis species diverged before their hosts, while the asynchronous timing of population declines suggests host shifts. Our results have uncovered complex patterns of genetic variation influenced by multiple factors that shaped the adaptation of Pneumocystis populations during their spread across mammals.IMPORTANCE Understanding how natural pathogen populations evolve and identifying the determinants of genetic variation are central issues in evolutionary biology. Pneumocystis, a fungal pathogen which infects mammals exclusively, provides opportunities to explore these issues. In humans, Pneumocystis can cause a life-threatening pneumonia in immunosuppressed individuals. In analysis of different Pneumocystis species infecting humans, rats, and mice, we found that there are high infection rates and that natural populations maintain a high level of genetic variation despite low levels of recombination. We found no evidence of population structuring by geography. Our comparisons of the times of divergence of these species to their respective hosts suggest that Pneumocystis may have undergone recent host shifts. The results demonstrate that Pneumocystis strains are widely disseminated geographically and provide a new understanding of the evolution of these pathogens.

Keywords: evolutionary biology; genetic diversity; genetic recombination; pneumonia; population structure.

FIG 1

Phylogenomic analysis reveals a lack of global population structuring in Pneumocystis jirovecii. Each symbol represents a single P. jirovecii isolate. The color-coding scheme indicates the region (city or state and country) in which the isolate was obtained. There is no clustering by location but instead an intermixing of isolates from diverse geographic locations, consistent with a lack of geographic structure. The phylogenetic tree was estimated from genome-wide SNP differences of Pneumocystis jirovecii isolates sequenced in this project (832,669 segregating sites).

FIG 2

Comparison of population genetic statistics, nucleotide diversity (π), Tajima’s D, Fay and Wu’s H, and linkage disequilibrium (r²) in three Pneumocystis species. (A) Box plots of genome-wide nucleotide diversity. The data show elevated median levels of nucleotide diversity in P. jirovecii relative to P. carinii and P. murina populations. (B and C) Box plots of summaries of site frequency spectra (Tajima's D [B] and Fay and Wu’s H [C]) data show that these statistics are skewed toward low-frequency variants in P. carinii and, to a lesser extent, in P. murina but not in P. jirovecii (negative Tajima’s D values), which suggests that these species experienced different demographic events. (D) Box plots of squared correlations (r²) between pairs of SNPs indicate strong signals of linkage disequilibrium (LD) in P. jirovecii and P. carinii. Signs above the bars indicate significant differences between species determined by the Mann-Whitney U test (***, P value of <2.2 × 10⁻¹⁶). The boxes indicate medians and interquartile ranges, whiskers indicate 95% values, and additional points in each box plot represent outliers.

FIG 3

Decay of linkage disequilibrium (LD) as a function of distance. LD levels measured by the square of the correlation coefficient between two markers (r²) were calculated for all pairs of biallelic SNPs within a 1-kb genomic window and averaged. LD decay levels were below 0.2 at 123 and 114 kb for Pneumocystis jirovecii and P. carinii, respectively (shown as gray vertical lines). The low rates of LD decay suggest low rates of recombination in Pneumocystis, which would be at least 2-fold lower than Saccharomyces cerevisiae recombination rates (32).

FIG 4

Estimation of evolutionary rates using P. jirovecii mitogenome data. (A) Midpoint rooted maximum likelihood phylogenetic tree of 15 full-length mitochondrial genomes, with bootstrap values shown on the branches. The tree tips indicate the collection date (year) and the sample identifier; the color code indicates the country of origin for the sample. (B) Root-to-tip distances (numbers of substitutions per site × 10⁻³) significantly correlate with the collection dates, indicating the presence of temporal signals. Calculations of regression of distances against the dates were performed using Murray’s R scripts (46). Statistical significance data are based on 1,000 random permutations.

FIG 5

Extended Bayesian skyline plot (EBSP) results for populations of two Pneumocystis species. The gray areas represent the upper and lower 95% highest posterior density (HPD) data, and the black dashed line represents the median estimate of 95% highest-posterior-density bounds. (A) Each graph represents the results of an EBSP analysis for Pneumocystis jirovecii (A1) and for P. jirovecii (A2), with the timeline restricted to 500,000 years before now. The purple dot indicates the approximate time of population decline. (B) Pneumocystis carinii (B1) and P. carinii (B2). The timeline was restricted to 20,000 years before now. The y axis shows the effective population size (Ne) per generation time (~4.5 days).

FIG 6

Overview of the timing and evolution of Pneumocystis and their mammalian hosts. (A) Phylogeny and coalescent divergence time estimates of Pneumocystis and their mammalian hosts were determined using nuclear genic alignments. The divergence time estimates (in millions of years) are displayed at the nodes. (B) Time line of events in Pneumocystis evolution relative to the mammal evolutionary history. For all Pneumocystis species, the time of divergence precedes the time of divergence of their respective hosts, which is inconsistent with the hypothesis that Pneumocystis coevolved with its hosts.