Original Research ARTICLE
Diversity of Pneumocystis jirovecii during Infection Revealed by Ultra-Deep Pyrosequencing
- 1Laboratoire de Parasitologie-Mycologie, Groupe Hospitalier Saint-Louis-Lariboisière-Fernand-Widal, Assistance Publique Hôpitaux de Paris, Hôpital Saint-Louis, Paris, France
- 2Université Paris Diderot, Sorbonne Paris Cité, Paris, France
- 3Unité de Mycologie Moléculaire, Département de Mycologie, Centre National de Référence Mycoses Invasives et Antifongiques, Institut Pasteur, Paris, France
- 4Centre National de la Recherche Scientifique CNRS URA3012, Paris, France
- 5Laboratoire de Microbiologie, Groupe Hospitalier Saint-Louis-Lariboisière-Fernand-Widal, Assistance Publique Hôpitaux de Paris, Hôpital Saint-Louis, Paris, France
Pneumocystis jirovecii is an uncultivable fungal pathogen responsible for Pneumocystis pneumonia (PCP) in immunocompromised patients, the physiopathology of which is only partially understood. The diversity of the Pneumocystis strains associated with acute infection has mainly been studied by Sanger sequencing techniques precluding any identification of rare genetic events (< 20% frequency). We used next-generation sequencing to detect minority variants causing infection, and analyzed the complexity of the genomes of infection-causing P. jirovecii. Ultra-deep pyrosequencing (UDPS) of PCR amplicons of two nuclear target region [internal transcribed spacer 2 (ITS2) and dihydrofolate reductase (DHFR)] and one mitochondrial DNA target region [the mitochondrial ribosomal RNA large subunit gene (mtLSU)] was performed on 31 samples from 25 patients. UDPS revealed that almost all patients (n = 23/25, 92%) were infected with mixtures of strains. An analysis of repeated samples from six patients showed that the proportion of each variant change significantly (by up to 30%) over time on treatment in three of these patients. A comparison of mitochondrial and nuclear UDPS data revealed heteroplasmy in P. jirovecii. The recognition site for the homing endonuclease I-SceI was recovered from the mtLSU gene, whereas its two conserved motifs of the enzyme were not. This suggests that heteroplasmy may result from recombination induced by unidentified homing endonucleases. This study sheds new light on the biology of P. jirovecii during infection. PCP results from infection not with a single microorganism, but with a complex mixture of different genotypes, the proportions of which change over time due to intricate selection and reinfection mechanisms that may differ between patients, treatments, and predisposing diseases.
Pneumocystis jirovecii is an ascomycete fungus that specifically infects humans (Cushion, 2010; Gigliotti and Wright, 2012) and causes Pneumocystis pneumonia (PCP), mostly in immunocompromised patients, including HIV-positive, solid organ transplant, and cancer/hematology patients, but also in adults and children with various other underlying conditions (Pagano et al., 2002; Roblot et al., 2003; Catherinot et al., 2010; Wissmann et al., 2010; Reid et al., 2011; Mori and Sugimoto, 2012; Tasaka and Tokuda, 2012). P. jirovecii thrives at the surface of alveolar pneumocytes, but is unable to grow on standard artificial media. It can, however, be amplified in vitro in an air-liquid interface culture system (Schildgen et al., 2014). These characteristics have made it difficult to study its genetic diversity, complexity, and evolution in humans. In addition, the genome sequence of P. jirovecii is only recently publicly available (Cissé et al., 2012; Cushion and Keely, 2013).
Since the 1990s, specific tools have been used to study the genetic diversity of P. jirovecii. In addition to cloning and sequencing (Tsolaki et al., 1996), single-strand conformation polymorphism (Hauser et al., 2001; Hauser, 2004), single-nucleotide extension assays (Esteves et al., 2011; Alanio et al., 2015), multi-locus sequencing typing (MLST) (Maitte et al., 2013), and short tandem repeat (STR)-based genotyping methods (Parobek et al., 2014; Gits-Muselli et al., 2015) have shown that mixtures of genotypes, differing in at nuclear or mitochondrial loci, are present during active infection. It has been estimated that up to 70% of patients with PCP harbor multiple genotypes, as assessed in various populations and by diverse techniques targeting nuclear or mitochondrial genes (Hauser et al., 2001; Parobek et al., 2014; Gits-Muselli et al., 2015). Studies based on Sanger sequencing have also reported the presence of mixed P. jirovecii sequences, but in a lower proportion of patients (about 30%) due to the well-known poor sensitivity of the Sanger method for detecting minority alleles (Angulo et al., 2010).
The genetic diversity of P. jirovecii has mostly been investigated in epidemiological studies, in which little attention was paid to the differences between nuclear and mitochondrial markers. The mitochondrial ribosomal RNA large subunit gene is currently used as a PCR target for PCP diagnosis (Meliani et al., 2003; Alanio et al., 2011, 2016; Botterel et al., 2012) and for typing based on MLST (Maitte et al., 2013). However, the nuclear and mitochondrial genomes can behave differently. Indeed, mitochondrial DNA can undergo more recombination events and acquire more mutations over time than nuclear DNA (Fritsch et al., 2014). The mitochondrial genome of P. jirovecii has recently been sequenced and compared to those of other Pneumocystis species (Ma et al., 2013). Although the mitochondrial genome of Pneumocystis carinii and Pneumocystis murina are thought to be linear, that of P. jirovecii appears to consist of a double-stranded circular mitochondrial DNA (mtDNA) molecule (Ma et al., 2013) that is autonomously replicated, transcribed, and mostly encodes proteins involved in mitochondrial respiration (Chen and Butow, 2005). In fungi, the mitochondrial genome mostly displays uniparental inheritance (Basse, 2010) and is autonomously replicated and transcribed (Chatre and Ricchetti, 2014).
We therefore investigated the diversity of P. jirovecii sequences, using ultra-deep pyrosequencing (UDPS) to detect minority variants and to compare the dynamics of the mitochondrial and nuclear genomes in the context of infection.
Materials and Methods
Saint-Louis Hospital, Paris, France, is a 650-bed tertiary university hospital with major clinical activities in hematology, renal transplantation, and oncology. This study was a non-interventional study with no change in the usual procedures. Biological material and clinical data were obtained only for standard diagnostic following physicians' prescriptions with no specific sampling. According to the French Health Public Law (CSP Art L1121-1.1), such protocol does not require approval of an ethics committee and is exempted from specific informed consent application.
Patients and Samples
We studied a collection of 31 respiratory samples from 25 patients with PCP treated at our university hospital (Saint Louis Hospital, Paris, France). Patients and samples were selected randomly from a collection of sputum or BAL fluids from patients with PCP and a high load of fungal DNA (Cq < 30, Alanio et al., 2011). Diagnosis was based on clinical (pulmonary infection in immunocompromised patients), radiological (interstitial pneumonia, ground glass opacities), and biological findings (ascus observed on immunofluorescence analyses of smears of respiratory material). DNA was extracted as previously described (Alanio et al., 2011). Our real-time quantitative PCR assay confirmed the high fungal load in all samples (>1900 EqTr/ml) (Alanio et al., 2011). The median age of the patients was 50 years (range: 1–80 years), and the sex ratio (M/F) was 5.25. The patients were HIV-positive (n = 12, 48%), had hematological malignancies (n = 6, 24%), were kidney transplant recipients (n = 4, 16%) or had another cause of immunodeficiency (n = 3, 12%). Two samples were collected from each of six patients (12 samples) with a median interval between samples of 7.5 (1–18) days. All samples had already been genotype with a recently described short-tandem repeat (STR) genotyping assay (Gits-Muselli et al., 2015). The six STR markers used were evenly distributed in the nuclear genome, as shown by alignment with the P. murina genome sequence (Gits-Muselli et al., 2015). The limit of detection for mixed STR markers was shown experimentally to be 2% (Gits-Muselli et al., 2015).
Selection of Target Genes and Amplicons
We were able to analyze amplicons of up to 320 bp in length in specific assays. The 5.8S-ITS2-28S (n = 120), DHFR (n = 60), and mtLSU (n = 52) sequence hits collected from Nucleotide Blast (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastSearch) were aligned, using Geneious 7.1 (Biomatters Ltd., New Zealand). We selected primers (Table S2) binding to conserved regions and allowing the amplification of sequences of about 300 bp in length, to maximize the number of polymorphic positions detected (Figures S1A–S3A). The primers were designed with Primer3web v 4.0.0 (http://primer3.ut.ee).
Specific Gene Amplification by UDPS
Each PCR primer consisted of a specific adaptor (5′-CGTATCGCCTCCCTCGCGCCATCAG-3′ for forward primers and 5′CTATGCGCCTTGCCAGCCCGCTCAG-3′ for reverse primers), a sample-specific multiplex identifier (MID) sequence, and a specific amplification sequence for PCR (Table S2).
We obtained three amplicons for each sample. Each amplicon (22.5 μl) was diluted with an equal volume of water and purified with AMPure beads (72 μl), before quantification with a Qubit® ds DNA HS Assay kit (Life Technologies) on a Qubit 3.0 instrument. For each gene, the PCR amplicons from 16 samples were pooled in equimolar concentrations and used to prepare an equimolar mixture of 107 molecules. This pool was subjected to clonal amplification on capture beads by emulsion PCR. In total, 500,000 beads enriched in DNA were deposited into the wells of a PicoTiter plate device and subjected to pyrosequencing in both directions with the GS Junior system (Roche). After a 10-h run, total reads were analyzed with GS amplicon variant analyzer (AVA) software (Roche), filtered and assigned, by demultiplexing, to the correct patient sample on the basis of MID correspondence.
Analysis of UDPS Data
AVA software (Roche) also aligned the generated sequence reads with the chosen reference sequence (JQ365725 for 5.8S-ITS2-28S, AF090368 for DHFR, JX499143 for mtLSU) and identified the proportion of polymorphism for each variable position in the sequences amplified. Insertions and deletions were not considered. Only already known SNPs were taken into account and haplotypes representing >0.5% of the sequences were selected for further analysis, to prevent the introduction of false polymorphic regions due to PCR errors (Brodin et al., 2013). Each individual combination of SNPs was selected and called haplotypes.
Reproducibility was determined by testing the same sample in two runs. The proportion of each of the variants sequenced varied by no more than 2%.
In all analyses, we studied one sample per patient (n = 25): the first sample collected. Specific analyses of repeated samples were performed separately.
Screening for the Presence of a Homing Endonuclease
We screened for the presence of conserved I-SceI motifs (GIGLILGDA and LAYWFMDDG) or motifs derived from this enzyme from Taphrina deformans (IVGLMLGDGY and LAVWISDDG) and for the recognition site of this enzyme (AAGTTACGCTAGGGATAACGGGGTAATATAGC) in the mitochondrial genome of P. jirovecii (JX499143), with Geneious v8.0 (Biomatters Ltd.).
Graphs and Statistical Analysis
Sequence analysis graphs were generated with Geneious v8.0 (Biomatters Ltd.). Hierarchical clustering was performed by Pearson's correlation analysis and average linkage clustering, with MeV v4.6.1 open-source genomic analysis software (The TM4 Development Group) obtained from www.tm4.org (Saeed et al., 2006). All graphs were plotted and statistical analyses (t- tests, One-way ANOVA, chi-square) were performed with Prism software v.6.0 (Graphpad).
UDPS Detected Numerous Variant Sequences
We used UDPS to analyze the diversity of three genetic sequences in 31 respiratory samples from 25 patients. We investigated two nuclear DNA target regions [partial 5.8S, internal transcribed spacer 2 and partial 28S gene (ITS2), GenBank accession n°JQ365725] and the [partial dihydrofolate reductase gene (DHFR), GenBank accession n°AF090368] and one mitochondrial DNA target region [partial mitochondrial ribosomal RNA large subunit gene (mtLSU), GenBank accession n°JX499143] (Figures S1A–S3A). The ITS2 locus was selected because this target DNA has already been used for genotyping P. jirovecii samples (de Boer et al., 2007; Le Gal et al., 2012) and is known to be present as a single copy in the P. jirovecii genome (Nahimana et al., 2000; Cissé et al., 2012; Cushion and Keely, 2013). The dihydrofolate reductase gene (DHFR, target gene of trimethoprim (diaminopyrimidine), one of the components of cotrimoxazole) was also selected because it was a single-copy gene and because mutations of this gene have already been described and potentially associated with the failure of cotrimoxazole prophylaxis (Nahimana et al., 2004; Queener et al., 2013). The mtLSU locus is the DNA target most widely used for diagnosis, it is known to be single-copy in the mitochondrial genome (Ma et al., 2013), with the mitochondria and mitochondrial genome being replicated several times within each organism. Two runs per gene were performed by multiplexing. After filtering, sequences were obtained and aligned with the corresponding reference sequence (with a mean of 6552 ± 3525, 8107 ± 2625, 4022 ± 1892 reads for ITS2, DHFR, and mtLSU, respectively).
Polymorphisms were observed at positions 135 (T135A), 230 (T230A), and 232 (G232A or G232C) of the 277 bp ITS2 amplicon, at positions 192 (T192C) and 218 (A218G) of the 300 bp DHFR amplicons, and at positions 84 (A84T or A84C) and 247 (C247T) of the 314-bp mtLSU amplicons (Figure 1, Table S1, Figures S1B–S3B). We used these polymorphisms to determine the haplotypes [combination of single nucleotide polymorphisms (SNPs)] of the 31 samples. Overall, seven variants (type of haplotypes) were found for ITS2 (TTG [WT], ATG, AAG, TAT, ATC, TTA, AAA), three for DHFR (TA [WT], CA, and TG) and six for mtLSU (AC [WT], TT, AT, CC, CT, TC) (Figure 1, Figures S1B–S3B).
Figure 1. Variants of the three Pneumocystis jirovecii target DNA regions observed in our collection of samples by ultra-deep pyrosequencing (UDPS). Analysis of all 31 samples led to the identification of seven variants for the internal transcribed spacer two region (ITS2), two for the partial dihydrofolate reductase gene (DHFR), and six for the mitochondrial ribosomal RNA large subunit gene (mtLSU), with respect to reference sequences (GenBank Accession number: JQ365725, AF090368, and JX499143, respectively). Of note, for ITS2, the numbering corresponds to the ITS2 region and not to that of the whole amplicon. For the ITS2 target, three polymorphic bases at positions 66, 161, and 163 of the ITS2 region resulted in the observation of seven variants: TTG, TTA, ATA, AAA, AAG, ATG, and ATC. For the DHFR gene, two polymorphic bases at positions 192 and 218 of the amplicon resulted in the observation of only three variants. For the mtLSU target, two polymorphic bases at positions 84 and 247 of the amplicon resulted in the observation of six variants: AG (WT), TT, AT, CC, CT, and TC. Dots indicate bases identical to the reference sequence. Polymorphic bases are shown in color. The size of the consensus base (sequence Logo) is proportional to the proportion of sequences polymorphic at the corresponding position. WT, wild-type sequence corresponding to the reference sequence; MT, mutated sequence corresponding to polymorphic variants.
We avoided bias in the interpretation of diversity, by analyzing the proportions of the variants in the 25 baseline samples, as repeat samples were available for six of the 25 patients. Overall, the number of different variants (Table 1, Figure 2) and the proportion of each variant (Figure 3) in each sample differed.
Table 1. Number of variants observed and resulting classifications of each sample for the three DNA target regions analyzed using UDPS and the six STR markers.
Figure 2. Distribution of the number of variants found in each sample for the mtLSU, ITS2, and DHFR DNA targets. Each of the 25 samples is represented as an independent dot (A). The number of variants was significantly different between the three genes (one-way ANOVA, p < 0.0001), with mtLSU and ITS having significantly more variants than DHFR. The bar represents the median for each target region. Each of the samples is represented by a line linking the number of variants found for each of the three genes (B). Similar patterns in the distribution of the variants can be seen for our 25 samples.
Figure 3. Proportion of each variant observed in 25 samples, by UDPS, for the mtLSU, ITS2, and DHFR DNA targets. The proportion of each variant after pyrosequencing analysis (AVA software, Roche) is depicted for each target DNA region, ITS2, DHFR, mtLSU, and for the 25 baseline samples. Each variant is represented by a different color and samples are grouped to allow the visualization of clusters based on the proportion of the variants (this grouping differs for the three targets). Pure samples (“p”) are shown for each target. Pure sequences (proportion of 100%) were observed in 5/25 samples for ITS2 (20%), 17/25 samples for DHFR (68%), and 2/25 samples for mtLSU (8%). The number of times each variant was recovered and the proportions of the variants differed between samples, with some clusters observed. Based on the combined data for the two nuclear DNA targets (ITS2, DHFR), pure sequences were observed in 5/25 (20%) samples (samples 090, 149, 175, 192, 252). Two of these samples were also identified as pure (“p”) for the mtLSU DNA target (samples 090 and 192), whereas the other three samples harbored mixtures of sequences for mtLSU, consistent with mitochondrial heteroplasmy (“h”). Variants accounting for < 0.5% of all sequences and variants with unknown polymorphisms were not considered in this study.
Variation of the Number and Proportion of Variants, Revealing Possible Heteroplasmy
The distribution of the number of variants observed differed significantly between mtLSU, ITS2 and DHFR (Kruskal-Wallis, p < 0.0001, Figure 2A, Figure S4A). Up to four, two and six variants/sample were observed for ITS2, DHFR, and mtLSU, respectively. No technical bias was observed, because there was no significant correlation (p > 0.05) between the number of variants and fungal load (based on quantification cycle or Cq value) or total number of reads sequenced (Figure S4B). Clusters of samples based on the number of variants for the three targets could be defined (Table 1, Figure 2B). These clusters were not significantly associated with the patient's background (p > 0.05).
The proportion of the total number of reads sequenced for each target corresponding to each of the variants was determined for each of the 25 samples (Figure 3). Samples were classified as pure (single variant on UDPS analysis) or mixed (≥2 variants) sequences. Pure sequences were observed in only 5/25 samples for ITS2 (20%), 17/25 samples for DHFR (68%), and 2/25 (8%) of samples for mtLSU (Table 1, Figure 3). Some of the samples could be clustered on the basis of the proportions of variants, but the pattern of clustering differed between target genes (Figure 4).
Figure 4. Heatmap of the proportions of the different variants in the 31 samples analyzed by UDPS, for the mtLSU, ITS2, and DHFR DNA targets. The 31 samples were clustered by hierarchical clustering analysis based on the proportions of the variants recovered. Samples from the six patients from whom repeat samples were obtained are indicated by a colored box (P14 in red, P2 in green, P24 in orange, P25 in black, and P3 in blue). The clustering analysis showed that, for P2, P18, and P3, the proportion of variants varied little between the two samples, which clustered together. By contrast, for P14, P24, and P25, significant differences in the proportions of variants were observed, leading to the two samples from the same patient being assigned to different clusters. Gray indicates an absence of recovery, blue indicates a low proportion (< 10%) and yellow indicates a high proportion (>70%), as indicated by the scale bar at the top of each heatmap.
Based on the combined UDPS data for the two nuclear target genes (ITS2, DHFR), pure sequences were observed in 5/25 (20%) patients (samples #090, 149, 175, 192, 252) (Figure 2). The genotyping results obtained with the six STR nuclear markers (Gits-Muselli et al., 2015) were consistent with purity (a single allele detected for each of the six markers) for three of these five samples (#090, 192, 252) (Table 1). Two of these five samples (#090, 192) harbored a single mtLSU variant. The other three (#252, 149, 175) had mixed mtLSU sequences, with three, six and two mtLSU variants, respectively (Table 1). These results (polymorphic mtLSU sequences with unique nuclear markers) were suggestive of heteroplasmy.
Homing Endonuclease May Be Present in Pneumocystis jirovecii
We searched for the presence of a homing endonuclease site in the mtLSU sequence, because such heteroplasmy could have been acquired through mtDNA recombination events initiated by intron homing (Haugen et al., 2005). Analysis of the complete mtLSU sequence (JX499143) revealed the presence of a recognition site 81.25% identical to that described for the I-SceI of Saccharomyces cerevisiae (Dujon, 1980; Jacquier and Dujon, 1985). Only six of the 32 nucleotides differed from those present in the motif described for I-SceI (Figure 5) (Belfort and Roberts, 1997). The two conserved motifs (GIGLILGDA and LAYWFMDDG) of I-SceI were not found in the mtLSU DNA sequence, even if we used the motifs present in T. deformans (IVGLMLGDGY and LAVWISDDG), which, like Pneumocystis, belongs to the Taphrinomycotina subphylum (Cissé et al., 2014; Almeida et al., 2015).
Figure 5. Alignment of the recognition site for the homing endonuclease I-SceI with the corresponding region in the P. jirovecii mtLSU gene. The I-SceI homing endonuclease has a specific recognition site for double-strand DNA break and intron introduction. A putative recognition site for this enzyme was identified between positions 14,635 and 14,666 (JX499143). Nucleotides identical to the consensus sequence are shown in gray, whereas polymorphic positions are shown in color. Solid- and dotted-line arrows indicate cleavages on the presented and complementary strands, respectively. The triangle indicates the intron insertion site (adapted from Belfort and Roberts, 1997).
Patient-Based Analysis of P. jirovecii Diversity
The characteristics of the 25 patients are summarized in Table 2. Sex ratio was 7.3 and median age was 54 years. The patients were HIV-positive (n = 15), patients with hematological malignancies (n = 7), renal transplant recipients (n = 4), and had other causes of immunosuppression (n = 5). The country of birth was recovered in 21 patients and was France (n = 14), Haiti (n = 2), North Africa (n = 2), Romania (n = 1), Vietnam (n = 1), Comoro Islands (n = 1). In total, no significant correlation (p > 0.05) was found between the presence of the different variants for each gene or for the combination of the three genes and the gender, the age, the background, the country of birth of the patients and the evolution at month 3.
For the six patients with repeated samples (median time between the first and second samples: 7.5 days, range: 1–18 days), the proportions of the variants varied by more than 2% (cutoff for reproducibility) in all patients, and by more than 10% in two patients (P14, P25) (Table 3). In addition to the differences in the proportions of the variants present, some of the variants were present in one sample, and absent from the other in five patients (P2, P3, P14, P24, and P25), although the proportion of this variant was below 2% when it was present in some cases (P2, P3, P14, P24). Even for samples recovered only a few days apart (2 days for P2, and 1 day for P14), the proportions of the variants varied by up to 24.9% between samples (ITS2 and TAT for P2, ITS2 and TTG for P14) (Table 3). In P25, one DHFR variant (CA) accounted for 30.4% of the variants present 8 days after the first sample was collected and cotrimoxazole treatment was initiated (Table 3). Clustering on the basis of the proportion of variants for the three targets revealed that the two samples from P2, P3, and P18 clustered together, whereas those from P14, P24, and P25 did not (Figure 4).
Table 3. Description of the proportion of variants in the repeated samples from six patients based on three DNA target regions analyzed using UDPS.
Our next-generation sequencing results indicate that Pneumocystis infections are mostly mixed infections. UDPS of nuclear gene targets revealed that 80% of the samples harbored mixed sequences, and this proportion increased to 92% (23/25) when a mitochondrial DNA target was added to the analysis.
Primary infection occurs early in life, with about 85% of children displaying anti-Pneumocystis antibody-mediated immune responses by the age of 18 months (Meuwissen et al., 1977; Pifer et al., 1978) and the detection of fungal DNA in non-invasive respiratory samples (Vargas et al., 2001). Immunocompetent or immunosuppressed individuals act as a reservoir for P. jirovecii, which infects immunocompromised patients via airborne transmission (Gigliotti et al., 2003; Chabé et al., 2004; Choukri et al., 2010), as demonstrated by epidemiological studies of numerous PCP outbreaks in a context of intrahospital transmission, particularly in renal transplant units (de Boer et al., 2011). The high predominance of mixed infection observed in our study may reflect continuous exposure throughout life (Gigliotti and Wright, 2012), together with the active multiplication of a subset of strains during immunosuppression. Some strains may be less capable of active multiplication than others in specific settings and in a given patient. This does not exclude the possibility that reactivation from dormancy contributes to the diversity observed during infection, because differences have been observed in the geographic distributions of P. jirovecii genotypes (Alanio et al., 2015). From an epidemiologic standpoint, these findings demonstrate that mitochondrial genes, including the mtLSU gene in particular, should be used with caution for the MLST genotyping of P. jirovecii because the mitochondrial and nuclear genomes may evolve differently.
Different mechanisms are thought to be involved in mutations of the mitochondrial genome. Mitochondrial DNA is subjected to more rounds of replication than nuclear DNA, and this may result in a higher frequency of point mutations. Mutations may also accumulate more easily due to higher levels of exposure to reactive oxygen species and, finally, the DNA repair machinery of the mitochondria has been shown to be less efficient than that of the nucleus (Song et al., 2005). Interestingly, five samples harbored pure nuclear (ITS and DHFR) sequences but mixed mitochondrial sequences, suggesting that heteroplasmy had occurred. Of course, we cannot completely rule out the possibility of polymorphism in other nuclear sequences. However, the analysis of six additional nuclear markers by STR typing (high level of discrimination and widespread distribution within the P. jirovecii genome; Gits-Muselli et al., 2015) was consistent with purity (unique genotype) for three samples. Heteroplasmy has already been reported in fungi, mostly for filamentous fungi (Barr et al., 2005). However, heteroplasmy seems to be rare in most fungi, due to uniparental mitochondrial transmission (Barr et al., 2005). In laboratory conditions, heteroplasmy is controlled and used to study mtDNA recombination (Basse, 2010). By contrast, heteroplasmy has rarely been studied in natural populations. Mitochondrial gene variation between isolates has been demonstrated by cloning in arbuscular mycorrhizal fungi (Rhizophagus irregularis) (Beaudet et al., 2015), with the identification of up to 20 variants per isolate for various mitochondrial target regions, all different from mtLSU. This aspect has been studied in natural populations of only two medically important fungi, Candida albicans (Anderson et al., 2001) and Cryptococcus gattii (Xu et al., 2009), to our knowledge. In both cases, the authors used PCR followed by Sanger sequencing and strain culture, assuming that the culture was clonal and therefore had only one nuclear genome. As expected, given the greater sensitivity of UDPS, we detected a higher degree of polymorphism within P. jirovecii mtLSU. Moreover, we had to control for the uniqueness of the nuclear genome in the absence of a simple culture system. Heteroplasmy was, therefore, expected in UDPS with a 0.5% cutoff value for variant selection. The uncultivable nature of P. jirovecii complicates studies of this phenomenon and the interpretation of its importance. However, recently developed culture methods may facilitate studies of the dynamics of heteroplasmy (Schildgen et al., 2014). However, homoplasmy is widely described in various organisms and controlled by genetic mechanisms (Christie et al., 2015). The observation of heteroplasmy in P. jirovecii could be the consequence of the lack of the control mechanisms that tend to decrease heteroplasmy rates are not present in P. jirovecii, or their ineffectiveness during the active phase of disease. P. jirovecii could also regularly or occasionally inherit mitochondria from both parents. In addition, a genetic phenomenon called intron homing responsible for recombination in the mitochondrial genome could lead to heteroplasmy.
Homing endonuclease genes (HEGs) are widespread in fungal genomes and are generally located in self-splicing introns (Basse, 2010). A recognition site for the homing endonuclease I-SceI was found in the S. cerevisiae mtLSU gene (Dujon, 1980; Jacquier and Dujon, 1985). This suggests that heteroplasmy in P. jirovecii may result from recombination due to the double-strand breaks created by such enzymes in mitochondrial DNA. However, the classical LAGLIDADG (Belfort and Roberts, 1997) motif and the conserved motifs GIGLILGDA and LAYWFMDDG (Chevalier et al., 2005) were not detected in the mtLSU gene, even if the mitochondrial homolog from T. deformans (the closest relative of Pneumocystis; Cissé et al., 2013, 2014; Almeida et al., 2015) was used as the reference (GenBank: CAUL01000001). This suggests that there may be as yet undiscovered homologs of I-SceI endonucleases in P. jirovecii.
Several variants, including those of the mtLSU target, were clearly present at the time of disease diagnosis, when fungal load was high. The observed diversity may also be due to mutations acquired during active growth of the fungus during immunosuppression and in the days or weeks before diagnosis. We cannot exclude the possibility of a peculiar variant having a selective advantage over another variant, resulting in changes in their relative proportions overtime. A stable final equilibrium in the proportion of variants may be reached in some patients, but not in others. Indeed, we observed variations in the proportions of variants in analyses of results for the serial samples obtained from some patients. This variation may also be due to differences in the sampling method, because P. jirovecii isolates recovered from BAL or sputum have already been shown to be different from those recovered from lung tissue (Helweg-Larsen et al., 2001). One patient (P25) with AIDS acquired a mutation of the DHFR target DNA (T192C) in about 30% of all DHFR sequences 8 days after the initiation of cotrimoxazole treatment (Table S1). This mutation corresponds to a T1158C synonymous substitution (called T312C in Nahimana et al., 2004) with no effect on protein function and not reflecting the acquisition of a resistance mutation. The other coding regions of this gene have not been studied and may contain other mutations. However, the observed variant accounted for only 30% of all variants, and complete recovery with no clinical resistance to cotrimoxazole was observed. The detection of this variant only 8 days after the introduction of cotrimoxazole treatment could suggest acquisition of a mutation during proliferation in the presence of the drug or selection and proliferation of this under treatment pressure.
Our study provides additional evidence for the complex natural history of P. jirovecii infections, for which the three Henle-Koch's postulates long used to evaluate the causal relationship between an infectious agent and a clinically defined disease can't be applied (Evans, 1976; Falkow, 1988). We show that these infections are mostly associated with mixtures of strains/genotypes that were probably acquired over the course of the patient's life. This finding is consistent with the current view of the natural history of these infections: reservoir in individual with various backgrounds (immunocompetent infants and adults or individuals with chronic pulmonary diseases such as chronic obstructive pulmonary disease or with relative immunosuppression such as children and pregnant women; Peterson and Cushion, 2005) and direct airborne transmission (Cushion, 2010; Gigliotti and Wright, 2012). However, we still know little about the subjects carrying P. jirovecii but not developing active infection (Hauser et al., 2000; Morris and Norris, 2012), as the low fungal load in respiratory samples from these individuals precludes such analyses, particularly for single-copy nuclear DNA targets.
Conceived and designed the experiments: AA, SB. Performed the experiments: MG, SM. Analyzed the data: AA, FD, SB. Contributed reagents/materials/analysis tools: SM, FD. Wrote the paper: AA, FD, SB. Major criticism of the manuscript: FD, SB.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank François Simon from the Microbiology Department of Saint Louis Hospital, Paris for allowing us to use the GS Junior device (Roche Diagnostics), François Collin (Roche Diagnostics) for technical support, and all the physicians involved in the collection of samples and patient care at Saint-Louis Hospital (Prof. Anne Bergeron, Prof. Marie-Noelle Péraldi, Prof. Elie Azoulay, Prof. Jean-Michel Molina, Dr. Nathalie de Castro, and Dr. Blandine Denis).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fmicb.2016.00733
Figure S1. Alignment of ITS2 sequences recovered from GenBank (n = 120) (A) and of the amplicons recovered from our patients, with the polymorphic bases highlighted (B).
Figure S2. Alignment of DHFR sequences recovered from GenBank (n = 60) (A) and of the amplicons recovered from our patients, with the polymorphic bases highlighted (B).
Figure S3. Alignment of mtLSU sequences recovered from GenBank (n = 52) (A) and of the amplicons recovered from our patients, with the polymorphic bases highlighted (B).
Figure S4. Distribution of the number of variants found in each sample for mtLSU, ITS2, and DHFR DNA targets, as determined by UDPS and with six nuclear STR markers (A). Distribution of fungal load (left panel) and the number of reads for each patient (one dot per patient, n = 31), classified by the maximum number of variants recovered per sample (1–4 for ITS, 1–2 for DHFR, and 1–6 for mtLSU) (B).
Table S1. Correspondence between the position of the observed polymorphic bases described in our samples and their position in the reference sequence, with the corresponding amino acid substitution.
Table S2. Primers and PCR conditions used.
Alanio, A., Desoubeaux, G., Sarfati, C., Hamane, S., Bergeron, A., Azoulay, E., et al. (2011). Real-time PCR assay-based strategy for differentiation between active Pneumocystis jirovecii pneumonia and colonization in immunocompromised patients. Clin. Microbiol. Infect. 17, 1531–1537. doi: 10.1111/j.1469-0691.2010.03400.x
Alanio, A., Hauser, P. M., Lagrou, K., Melchers, W. J., Helweg-Larsen, J., Matos, O. et al. (2016). ECIL guidelines for the diagnosis of Pneumocystis jirovecii pneumonia in patients with haematological malignancies and stem cell transplant recipients. J. Antimicrob. Chemother. (in press).
Alanio, A., Olivi, M., Cabaret, O., Foulet, F., Bellanger, A.-P., Millon, L., et al. (2015). Correlation between Pneumocystis jirovecii mitochondrial genotypes and high and low fungal loads assessed by single nucleotide primer extension assay and quantitative real-time PCR. J. Eukaryot. Microbiol. 62, 650–656. doi: 10.1111/jeu.12222
Almeida, J. M. G. C. F., Cissé, O. H., Fonseca, Á., Pagni, M., and Hauser, P. M. (2015). Comparative genomics suggests primary homothallism of Pneumocystis species. mBio 6, e02250–e02214. doi: 10.1128/mBio.02250-14
Anderson, J. B., Wickens, C., Khan, M., Cowen, L. E., Federspiel, N., Jones, T., et al. (2001). Infrequent genetic exchange and recombination in the mitochondrial genome of Candida albicans. J. Bacteriol. 183, 865–872. doi: 10.1128/JB.183.3.865-872.2001
Angulo, B., García-García, E., Martínez, R., Suárez-Gauthier, A., Conde, E., Hidalgo, M., et al. (2010). A commercial real-time PCR kit provides greater sensitivity than direct sequencing to detect KRAS mutations: a morphology-based approach in colorectal carcinoma. J. Mol. Diagn. 12, 292–299. doi: 10.2353/jmoldx.2010.090139
Beaudet, D., de la Providencia, I. E., Labridy, M., Roy-Bolduc, A., Daubois, L., and Hijri, M. (2015). Intraisolate mitochondrial genetic polymorphism and gene variants coexpression in arbuscular mycorrhizal fungi. Genome Biol. Evol. 7, 218–227. doi: 10.1093/gbe/evu275
Botterel, F., Cabaret, O., Foulet, F., Cordonnier, C., Costa, J.-M., and Bretagne, S. (2012). Clinical significance of quantifying Pneumocystis jirovecii DNA by using real-time PCR in bronchoalveolar lavage fluid from immunocompromised patients. J. Clin. Microbiol. 50, 227–231. doi: 10.1128/JCM.06036-11
Brodin, J., Mild, M., Hedskog, C., Sherwood, E., Leitner, T., Andersson, B., et al. (2013). PCR-induced transitions are the major source of error in cleaned ultra-deep pyrosequencing data. PLoS ONE 8:e70388. doi: 10.1371/journal.pone.0070388
Catherinot, E., Lanternier, F., Bougnoux, M.-E., Lecuit, M., Couderc, L.-J., and Lortholary, O. (2010). Pneumocystis jirovecii Pneumonia. Infect. Dis. Clin. North Am. 24, 107–138. doi: 10.1016/j.idc.2009.10.010
Chabé, M., Dei-Cas, E., Creusy, C., Fleurisse, L., Respaldiza, N., Camus, D., et al. (2004). Immunocompetent hosts as a reservoir of pneumocystis organisms: histological and rt-PCR data demonstrate active replication. Eur. J. Clin. Microbiol. Infect. Dis. 23, 89–97. doi: 10.1007/s10096-003-1092-2
Chevalier, B., Monnat, R. J. Jr., and Stoddard, B. L. (2005). The LAGLIDADG homing endonuclease family, in Homing Endonucleases and Inteins Nucleic Acids and Molecular Biology, eds M. Belfort, D. W. Wood, B. L. Stoddard, and V. Derbyshire (Heidelberg: Springer Berlin Heidelberg), 33–47.
Choukri, F., Menotti, J., Sarfati, C., Lucet, J. C., Nevez, G., Garin, Y. J. F., et al. (2010). Quantification and Spread of Pneumocystis jirovecii in the surrounding air of patients with Pneumocystis Pneumonia. Clin. Infect. Dis. 51, 259–265. doi: 10.1086/653933
Christie, J. R., Schaerf, T. M., and Beekman, M. (2015). Selection against heteroplasmy explains the evolution of uniparental inheritance of mitochondria. PLoS Genet. 11:e1005112. doi: 10.1371/journal.pgen.1005112
Cissé, O. H., Almeida, J. M. G. C. F., Fonseca, Á., Kumar, A. A., Salojärvi, J., Overmyer, K., et al. (2013). Genome sequencing of the plant pathogen Taphrina deformans, the causal agent of peach leaf curl. mBio 4, e00055-e00013. doi: 10.1128/mbio.00055-13
Cissé, O. H., Pagni, M., and Hauser, P. M. (2012). De novo assembly of the Pneumocystis jirovecii genome from a single bronchoalveolar lavage fluid specimen from a patient. mBio 4, e00428-e00412. doi: 10.1128/mBio.00428-12
Cissé, O. H., Pagni, M., and Hauser, P. M. (2014). Comparative genomics suggests that the human pathogenic fungus Pneumocystis jirovecii acquired obligate biotrophy through gene loss. Genome Biol. Evol. 6, 1938–1948. doi: 10.1093/gbe/evu155
Cushion, M. T. (2010). Are members of the fungal genus Pneumocystis (a) Commensals, (b) Opportunists, (c) Pathogens, or (d) All of the Above? PLoS Pathog. 6:e1001009. doi: 10.1371/journal.ppat.1001009
de Boer, M. G. J., Bruijnesteijn van Coppenraet, L. E. S., Gaasbeek, A., Berger, S. P., Gelinck, L. B. S., van Houwelingen, H. C., et al. (2007). An outbreak of Pneumocystis jiroveci pneumonia with 1 predominant genotype among renal transplant recipients: interhuman transmission or a common environmental source? Clin. Infect. Dis. 44, 1143–1149. doi: 10.1086/513198
de Boer, M. G. J., de Fijter, J. W., and Kroon, F. P. (2011). Outbreaks and clustering of Pneumocystis pneumonia in kidney transplant recipients: a systematic review. Med. Mycol. 49, 673–680. doi: 10.3109/13693786.2011.571294
Dujon, B. (1980). Sequence of the intron and flanking exons of the mitochondrial 21S rRNA gene of yeast strains having different alleles at the omega and rib-1 loci. Cell 20, 185–197. doi: 10.1016/0092-8674(80)90246-9
Esteves, F., Gaspar, J., De Sousa, B., Antunes, F., Mansinho, K., and Matos, O. (2011). Clinical relevance of multiple single-nucleotide polymorphisms in Pneumocystis jirovecii Pneumonia: development of a multiplex PCR-single-base-extension methodology. J. Clin. Microbiol. 49, 1810–1815. doi: 10.1128/JCM.02303-10
Gigliotti, F., Harmsen, A. G., and Wright, T. W. (2003). Characterization of transmission of Pneumocystis carinii f. sp. muris through immunocompetent BALB/c mice. Infect. Immun. 71, 3852–3856. doi: 10.1128/IAI.71.7.3852-3856.2003
Gits-Muselli, M., Peraldi, M.-N., de Castro, N., Delcey, V., Menotti, J., Guigue, N., et al. (2015). New short tandem repeat-based molecular typing method for Pneumocystis jirovecii reveals intrahospital transmission between patients from different wards. PLoS ONE 10:e0125763. doi: 10.1371/journal.pone.0125763
Hauser, P. M., Blanc, D. S., Bille, J., Nahimana, A., and Francioli, P. (2000). Carriage of Pneumocystis carinii by immunosuppressed patients and molecular typing of the organisms. AIDS 14, 461–463. doi: 10.1097/00002030-200003100-00022
Hauser, P. M., Blanc, D. S., Sudre, P., Senggen Manoloff, E., Nahimana, A., Bille, J., et al. (2001). Genetic diversity of Pneumocystis carinii in HIV-positive and -negative patients as revealed by PCR-SSCP typing. AIDS 15, 461–466. doi: 10.1097/00002030-200103090-00004
Helweg-Larsen, J., Lundgren, B., and Lundgren, J. D. (2001). Heterogeneity and compartmentalization of Pneumocystis carinii f. sp. hominis genotypes in autopsy lungs. J. Clin. Microbiol. 39, 3789–3792. doi: 10.1128/JCM.39.10.3789-3792.2001
Jacquier, A., and Dujon, B. (1985). An intron-encoded protein is active in a gene conversion process that spreads an intron into a mitochondrial gene. Cell 41, 383–394. doi: 10.1016/S0092-8674(85)80011-8
Le Gal, S., Damiani, C., Rouillé, A., Grall, A., Tréguer, L., Virmaux, M., et al. (2012). A cluster of Pneumocystis infections among renal transplant recipients: molecular evidence of colonized patients as potential infectious sources of Pneumocystis jirovecii. Clin. Infect. Dis. 54, e62–e71. doi: 10.1093/cid/cir996
Ma, L., Huang, D. W., Cuomo, C. A., Sykes, S., Fantoni, G., Das, B., et al. (2013). Sequencing and characterization of the complete mitochondrial genomes of three Pneumocystis species provide new insights into divergence between human and rodent Pneumocystis. FASEB J. 27, 1962–1972. doi: 10.1096/fj.12-224444
Maitte, C., Leterrier, M., Le Pape, P., Miegeville, M., and Morio, F. (2013). Multilocus sequence typing of Pneumocystis jirovecii from clinical samples: how many and which loci should be used? J. Clin. Microbiol. 51, 2843–2849. doi: 10.1128/JCM.01073-13
Meliani, L., Develoux, M., Marteau-Miltgen, M., Magne, D., Barbu, V., Poirot, J.-L., et al. (2003). Real time quantitative PCR assay for Pneumocystis jirovecii detection. J. Eukaryot. Microbiol. 50(Suppl.), 651. doi: 10.1111/j.1550-7408.2003.tb00670.x
Meuwissen, J. H., Tauber, I., Leeuwenberg, A. D., Beckers, P. J., and Sieben, M. (1977). Parasitologic and serologic observations of infection with Pneumocystis in humans. J. Infect. Dis. 136, 43–49. doi: 10.1093/infdis/136.1.43
Nahimana, A., Francioli, P., Blanc, D. S., Bille, J., Wakefield, A. E., and Hauser, P. M. (2000). Determination of the copy number of the nuclear rDNA and beta-tubulin genes of Pneumocystis carinii f. sp. hominis using PCR multicompetitors. J. Eukaryot. Microbiol. 47, 368–372. doi: 10.1111/j.1550-7408.2000.tb00062.x
Nahimana, A., Rabodonirina, M., Bille, J., Francioli, P., and Hauser, P. M. (2004). Mutations of Pneumocystis jirovecii dihydrofolate reductase associated with failure of prophylaxis. Antimicrob. Agents Chemother. 48, 4301–4305. doi: 10.1128/AAC.48.11.4301-4305.2004
Pagano, L., Fianchi, L., Mele, L., Girmenia, C., Offidani, M., Ricci, P., et al. (2002). Pneumocystis carinii pneumonia in patients with malignant haematological diseases: 10 years' experience of infection in GIMEMA centres. Br. J. Haematol. 117, 379–386. doi: 10.1046/j.1365-2141.2002.03419.x
Parobek, C. M., Jiang, L. Y., Patel, J. C., Alvarez-Martinez, M. J., Miro, J. M., Worodria, W., et al. (2014). Multilocus microsatellite genotyping array for investigation of genetic epidemiology of Pneumocystis jirovecii. J. Clin. Microbiol. 52, 1391–1399. doi: 10.1128/JCM.02531-13
Queener, S. F., Cody, V., Pace, J., Torkelson, P., and Gangjee, A. (2013). Trimethoprim resistance of dihydrofolate reductase variants from clinical isolates of Pneumocystis jirovecii. Antimicrob. Agents Chemother. 57, 4990–4998. doi: 10.1128/AAC.01161-13
Reid, A. B., Chen, S. C. A., and Worth, L. J. (2011). Pneumocystis jirovecii pneumonia in non-HIV-infected patients: new risks and diagnostic tools. Curr. Opin. Infect. Dis. 24, 534–544. doi: 10.1097/QCO.0b013e32834cac17
Roblot, F., Le Moal, G., Godet, C., Hutin, P., Texereau, M., Boyer, E., et al. (2003). Pneumocystis carinii pneumonia in patients with hematologic malignancies: a descriptive study. J. Infect. 47, 19–27. doi: 10.1016/S0163-4453(03)00038-0
Schildgen, V., Mai, S., Khalfaoui, S., Lüsebrink, J., Pieper, M., Tillmann, R. L., et al. (2014). Pneumocystis jirovecii can be productively cultured in differentiated CuFi-8 airway cells. mBio 5, e01186-14. doi: 10.1128/mbio.01186-14
Song, S., Pursell, Z. F., Copeland, W. C., Longley, M. J., Kunkel, T. A., and Mathews, C. K. (2005). DNA precursor asymmetries in mammalian tissue mitochondria and possible contribution to mutagenesis through reduced replication fidelity. Proc. Natl. Acad. Sci. U.S.A. 102, 4990–4995. doi: 10.1073/pnas.0500253102
Tasaka, S., and Tokuda, H. (2012). Pneumocystis jirovecii pneumonia in non-HIV-infected patients in the era of novel immunosuppressive therapies. J. Infect. Chemother. 18, 793–806. doi: 10.1007/s10156-012-0453-0
Tsolaki, A. G., Miller, R. F., Underwood, A. P., Banerji, S., and Wakefield, A. E. (1996). Genetic diversity at the internal transcribed spacer regions of the rRNA operon among isolates of Pneumocystis carinii from AIDS patients with recurrent pneumonia. J. Infect. Dis. 174, 141–156. doi: 10.1093/infdis/174.1.141
Vargas, S. L., Hughes, W. T., Santolaya, M. E., Ulloa, A. V., Ponce, C. A., Cabrera, C. E., et al. (2001). Search for primary infection by Pneumocystis carinii in a cohort of normal, healthy infants. Clin. Infect. Dis. 32, 855–861. doi: 10.1086/319340
Wissmann, G., Morilla, R., Martín-Garrido, I., Friaza, V., Respaldiza, N., Povedano, J., et al. (2010). Pneumocystis jirovecii colonization in patients treated with infliximab. Eur. J. Clin. Invest. 41, 343–348. doi: 10.1111/j.1365-2362.2010.02415.x
Keywords: ultra-deep pyrosequencing, Pneumocystis jirovecii, mixed infection, mitochondria, heteroplasmy, recombination
Citation: Alanio A, Gits-Muselli M, Mercier-Delarue S, Dromer F and Bretagne S (2016) Diversity of Pneumocystis jirovecii during Infection Revealed by Ultra-Deep Pyrosequencing. Front. Microbiol. 7:733. doi: 10.3389/fmicb.2016.00733
Received: 01 February 2016; Accepted: 02 May 2016;
Published: 24 May 2016.
Edited by:Frederic Lamoth, Lausanne University Hospital, Switzerland
Reviewed by:Anthony George Tsolaki, Brunel University London, UK
Alexey Porollo, Cincinnati Children's Hospital Medical Center, USA
Copyright © 2016 Alanio, Gits-Muselli, Mercier-Delarue, Dromer and Bretagne. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Alexandre Alanio, email@example.com