Genomic reconstruction of an azole-resistant Candida parapsilosis outbreak and the creation of a multilocus sequence typing scheme

Clinical case description and outbreak reconstruction

We report a prolonged outbreak of Candida parapsilosis that occurred across multiple healthcare institutions across Berlin, Germany during the period of 2018 - 2022 (Figure 1a). Cases were initially detected by physicians at Hospital #1 (name anonymized for privacy), who recognized an increase in isolates with a distinct antifungal susceptibility profile in their facility. The isolates were resistant to fluconazole (FLC^R) and voriconazole (VOR^R), but susceptible to other triazoles. After a period of approximately five months, during which no FLC^R and VOR^R infections were reported, infection isolates with this susceptibility profile were once again observed at Hospital #1. The outbreak was officially declared at that time (October 2019), and was formally investigated including environmental surveillance. As C. parapsilosis is known to colonize the ear, 21 ICU patients at Hospital #1 were screened with swabs of both ears, but none were positive for C. parapsilosis. To investigate possible transfer by staff, sixteen bedside stethoscopes were screened. Only one sample was positive, but the resulting isolate of C. parapsilosis was azole susceptible. In addition, opened and sealed hand creams used by staff were screened, but the pathogen was not detected. Hand sampling of staff was not performed as this is discouraged in Germany due to stigmatization.

• Download figure
• Open in new tab

Figure 1. Timeline and hospital network dynamics of a Candida parapsilosis outbreak in Berlin, Germany.

(a) Timeline histogram of the 33 outbreak cases. Dates are taken from the date of first positive culture. Color indicates the inpatient hospital at the time of isolation. A single case was identified in an ambulatory patient, which is indicated with “NA” for the inpatient hospital. (b) Inpatient stays of the outbreak cases represented as a network, where nodes (large circles) represent hospitals in the Berlin area and edges (lines) represent direct or indirect patient transfers between two institutions. Edge thickness depicts the relative number of transfers between two nodes. Individual patient case numbers are indicated within the node and along the edge. For the network, admission and discharge data beginning 12 months prior to the first detected outbreak and extending until the end of 2022 was utilized.

To prevent further cases, emphasis was placed on hygiene measures and particularly hand disinfection. Following implementation of the infection control measures, there was a nine month period without FLC^R and VOR^R strains reported. However, additional cases were detected in September 2021 and the start-stop pattern continued through 2022 (Fig 1a). Consultations with other hospitals and clinical laboratories in Berlin as well as the National Reference Center for Invasive Fungal Infections (NRZMyk) identified C. parapsilosis infections with this same antifungal susceptibility profile at five other institutions around Berlin. Subsequent epidemiological and genomic investigations revealed a total of 33 outbreak cases that are analyzed in the rest of this manuscript. The majority of outbreak cases were first reported as wound infections (39%), followed by infections associated with central venous catheters (29%) (Table 1). Outbreak isolates were all resistant to fluconazole with a median minimum in inhibitory concentration (MIC) of 64 mg/l and voriconazole (median MIC 1 mg/l), but susceptible to itraconazole (median MIC 0.125 mg/l) and posaconazole (median MIC 0.03 mg/l) (Table 1).

To understand the transmission dynamics of the outbreak, we reconstructed the temporal and spatial relationship of the outbreak using admission and discharge records. The reconstruction of inpatient admissions and transfers between medical centers in Berlin revealed that even though the largest fraction of outbreak cases (n = 20) were diagnosed at Hospital #1, the majority of these patients had been directly transferred to this facility from a major surgical facility (indicated as Hospital #0 in Figure 1b). Supporting nosocomial transmission, the majority of outbreak cases (31/33; 94%) had verifiable, temporal overlap with other cases, i.e. inpatient stays at a given hospital whose dates overlapped with >one other outbreak patient (Figure 1b). Despite lacking specific ward information for all outbreak cases, seven patient clusters involving a total of 22 patients had demonstrated spatial and temporal overlap, with outbreak patients confirmed to have been on the same ward at the same time as other outbreak cases, strongly supporting patient-to-patient spread (Figure 1b). In summary, we report a multiyear, sporadic outbreak of azole-resistant C. parapsilosis in Germany that was not limited to any single hospital but involved a network of hospitals in Berlin.

Genomic epidemiology of C. parapsilosis confirms the clonal nature of Berlin outbreak strains and reveals the presence of the outbreak strain on additional continents

To define the genomic epidemiology of C. parapsilosis in Germany and examine the genomic relationship between the putative outbreak cases from the hospitals in Berlin, we performed whole genome sequencing (WGS) on 68 isolates (33 outbreak isolates; 33 non-outbreak isolates) collected by the German National Reference Centre for Invasive Fungal Infections (NRZMyk) between 2016 and 2022. Sequenced isolates originated from 29 medical centers in 16 cities. 64% (n = 44) of the sequenced samples originated from Berlin, including samples suspected to be part of a local outbreak and 14 non-linked isolates from different medical centers. Publicly available sequence data for an additional 190 samples was also included to facilitate comparison of the German sequenced isolates with the global C. parapsilosis population, bringing the total number of samples in the dataset to 258 genomes. WGS samples analyzed had a mean depth of coverage of 146x (range 10-248x) after data quality control and mapping (Supplemental File 1).

When built into a whole genome phylogeny, samples from Germany were distributed throughout the tree, with no phylogenetic group being unique to a specific country or continent (Figure 2). The close genetic relationship between the outbreak strains was evident by their tight grouping and short branch lengths in the phylogeny, supporting the nosocomial transmission that was inferred from the spatial and temporal overlap in cases (Figure 1b). The quasi-clonal relationship between outbreak isolates was further confirmed by pairwise single nucleotide variant (SNV) analysis, where isolates in the Berlin outbreak clade differed by an average of 36 SNVs (Figure 3). In contrast, the average number of SNVs between non-outbreak isolates from other locations in Germany was significantly higher at 2,079 SNVs (p < 0.001 by student’s t-test). Further, the number of SNVs observed between duplicate sequencing of the same freezer stock of an isolate strain was 15 SNVs, indicating that the number of pairwise SNV differences between repeat sequencing of the same sample was similar to the genetic differentiation observed between different outbreak isolates. Strikingly, the public sequence data revealed samples from Turkey, South Korea, and Kuwait that were closely related to the outbreak strains from Berlin (Figure 2). These isolates were also reported to be fluconazole resistant ^27–29, highlighting the global spread of azole-resistant C. parapsilosis. These tightly related global strains differed from the Berlin outbreak strains by an average of 96 SNVs (range 18-191), compared to the mean 36 SNVs that the Berlin cluster differed from each other on average. Altogether, whole genome sequencing confirmed the close genetic relationship of the outbreak strains and its presence of this lineage in multiple other countries.

• Download figure
• Open in new tab

Figure 2. Whole genome phylogeny of Candida parapsilosis isolates from a global distribution and their azole susceptibility profiles and ERG11 polymorphisms.

Maximum-likelihood phylogeny constructed from SNV data from 38,567 variable positions. Tree is shown rooted at the midpoint. To increase tree readability, samples from the same sequencing study that showed little genetic differentiation were randomly downsampled to a maximum of 4 samples per sequencing project - with the exception of the samples from the samples described for the first time in our study. Branches are supported by ultrafast bootstrap values of >0.95 unless otherwise indicated. Where available, azole susceptibility data relative to EUCAST breakpoints and ERG11 polymorphisms are indicated. Strains that were resistant to one or more azole but where no polymorphisms in ERG11 were detected are marked with an asterisk.

• Download figure
• Open in new tab

Figure 3. SNV matrix of the German isolates sequenced in this study.

Heatmap showing the pairwise SNV differences between the novel German isolates described in this study. Metadata bands indicate (from top to bottom): fluconazole susceptibility, year of sample isolation, and whether the sample is from the city of Berlin (including both outbreak and non-outbreak strains) or another city in Germany.

ERG11 Y132F is the primary azole resistance mechanism among C. parapsilosis strains from Germany

Of the 68 German isolates newly described and sequenced in this study, 65% were resistant to one or more triazole (n = 44) (Supplemental File 1). ERG11 Y132F was the dominant resistant mechanism, with 36 resistant isolates possessing this mutation (n = 31 outbreak; n = 5 non-outbreak). Interestingly, the non-outbreak Y132F strains were in two distantly related groups of samples based on the phylogeny (Figure 2; Supplemental Figure 1). Each of these Y132F-containing phylogenetic groups contained isolates from geographically distant regions of Germany, suggesting that there are at least 3 genomically-distinct lineages of Y132F C. parapsilosis circulating widely throughout the country.

We also observed three closely related but unlinked samples from different cities in Germany that were pan-azole resistant and carried a G458S substitution in ERG11 (Figure 2; Supplemental Figure 1) that has been linked to clinical pan-azole resistance previously ^12,18,30,31. Also included in our dataset were eight isolates that were resistant to one more azoles, but did not contain any resistance-associated polymorphisms in ERG11. However, these isolates contained sequence polymorphisms in several genes that have been linked to drug resistance, including the efflux pumps MDR1 and CDR1B, the multidrug resistance gene regulator MRR1, and the transcriptional activator TAC1 (Supplemental Table 1). While a casual role for most of these polymorphisms in drug resistance have not been experimentally established, one resistant isolate harbored a L518F substitution in TAC1 experimentally proven to confer azole resistance ⁸. Among our newly described isolates and the public genomes we observed additional amino acid polymorphisms in ERG11 present in susceptible strains, including M178T, N283Y, and both resistant and susceptible strains R398I ^32–34. Together, we find that ERG11 Y132F is the dominant resistant mechanism among both outbreak and non-outbreak fluconazole-resistant strains in Germany and, surprisingly, discover that our outbreak strains are closely related to resistant samples previously described from East and West Asia.

Identification of markers for multilocus typing in C. parapsilosis

As outbreaks of C. parapsilosis are being increasingly reported, there is an urgent need for a rapid and reproducible tool capable of differentiating isolates. Motivated by this, we set out to identify a minimal set of loci that could be used for Sanger sequence-based typing of both general epidemiology and outbreak detection. Based on the whole genome phylogeny, we identified 14 distinct genetic clusters of C. parapsilosis in our dataset of 258 genomes (Figure 4a, b). The clusters contained an average of 17 samples (range 2-105) and were numbered to be concordant with the five major clades of C. parapsilosis described in ³⁵.

• Download figure
• Open in new tab

Figure 4. Comparison of whole genome and four loci phylogenies of C.parapsilosis.

(a-b) Maximum-likelihood whole genome SNV phylogeny with tips colored by country of origin (a) or genetic cluster based on whole-genome data (b). (c) Maximum-likelihood phylogeny constructed from ∼700 bp of sequence from CPAR2_101470, CPAR2_108720, CPAR2_212940, and CPAR2_808110. Tip color indicates the genetic cluster samples were assigned based on the whole genome phylogeny in panel b.

To identify loci with high discriminatory power, we performed a genome-wide screen of genetic variant profiles first using full-length gene coding sequences, followed by 750 bp regions of CDS with the high discriminatory power, and finally assessing these loci in combination. Using this approach, we identified four 750 bp loci in CPAR2_101470, CPAR2_108720, CPAR2_212940, and CPAR2_808110 (Table 2) that were capable of differentiating isolates into 31 unique sequence types with a Hunter’s Discriminatory Power of 0·87. As C. parapsilosis is a diploid organism, an appreciable fraction of the discriminatory power comes from accounting for whether observed sequence polymorphisms were heterozygous or homozygous. When nucleotide zygosity was not considered, we were only able to identify 22 sequence types and the discriminatory power was reduced to 0·81. The four selected loci showed a high degree of agreement with the whole genome phylogeny. Our typing scheme was able to correctly identify samples to their whole genome-defined genetic cluster with 97·6% sensitivity and 99·9% specificity and a phylogeny created from loci sequence largely replicated the structure of the whole genome tree (Figure 4c). To confirm that the selected loci were present as a single copy and the discriminatory power was not being inflated by the presence of multiple copies in the genome, we assessed them for genome-wide copy number changes. Using three different copy number variation softwares, we did not detect any CNV changes in these loci, either additional copies or loss of the region (Supplemental Figure 2). Based on orthology with Candida albicans and Saccharomyces cerevisiae, the identified loci have putative functions in the cytoskeleton, cell wall organization, and ion transport (Table 2).

Laboratory validation of in silico-identified markers

To validate our in silico typing, we designed primers to PCR amplify the four identified loci (Table 3). Using a representative 14 samples from our dataset, including four outbreak strains expected to have identical sequence type profiles, we amplified and Sanger sequenced the four loci. In each sample, we detected the expected variants from the in silico analysis by Sanger sequencing and did not discover any variants that were not in the whole genome data (Supplemental Table 2). Resultantly, we could correctly assign samples to the nine unique sequence types predicted from the in silico typing and the six genetic clusters assigned to them by the whole genome phylogeny. As expected, the four outbreak samples had identical sequences at all four loci and were assigned to the same sequence type. An important note is that Sanger sequencing providers may not report dual peaks caused by heterozygous sites using IUPAC nucleotide codes. As this was the case for us, this was resolved by manually recoding positions in the chromatograms where equal sized dual peaks were recorded (for example, a double peak reporting an A and a G were recoded to the corresponding IUPAC base of R). To facilitate the typing of novel isolates by future researchers and clinical labs, FASTA alignments containing nucleotide sequence from each sequence type of each of the four makers are provided in Supplemental File 2.

Demonstration of the utility of the typing scheme for rapidly identifying outbreak clones

To demonstrate the value of our typing scheme for the rapid identification of clones within an outbreak setting, we tested our typing loci on a set of seven novel isolates that were not part of the dataset used to create the typing scheme. These isolates were collected from three hospitals in Berlin and were suspected to be additional outbreak cases based on their antifungal susceptibility profile (FLC^R, VOR^R, POS^S, ITR^S) and their isolation from hospitals with previous outbreak cases. PCR amplification and Sanger sequencing of the four typing loci revealed that their sequence type was identical to the outbreak isolates analyzed by WGS, and strongly suggesting that they were additional outbreak cases (Supplemental Table 2). Taken together, we selected and validated four loci that can be used for sequence-based typing of C. parapsilosis and demonstrated their value in rapidly identifying closely related samples.

Outbreak cases and sample collection

The 68 novel strains of C. parapsilosis sequenced in this study were sent to the German National Reference Centre for Invasive Fungal Infections (NRZMyk) by German healthcare facilities and laboratories for species confirmation and susceptibility testing between 2014 and 2022 (Table 1). For tracking outbreak dynamics, admission and discharge records were obtained for all outbreak patients beginning January 2018, corresponding to one year before the first outbreak case, and extending to December 2022. Patient information was obtained from the sample submission form of the NRZMyk. Samples were de-identified via the assignment of an NRZ ID (e.g. NRZ-XXXX) which cannot be linked to a patient, hospital, or region. The Ethics Committee of the University Hospital Jena approved this study under registration number 2024-3255-Daten.

Strains submitted to the NRZMyk from outside Berlin were selected for whole genome sequencing (WGS) if they displayed resistance to one or more medical azole. For cities with azole resistant samples selected for sequencing, an additional number of azole susceptible strains were included from this city. In total 44 azole resistant strains (FLC^R, VOR^R, POS^S, or ITR^R) and 25 fully azole susceptible strains were analyzed. Strains were identified as C. parapsilosis prior to WGS by Sanger sequencing of the ITS locus. This locus was amplified with the forward primer V9G (5′- TTACGTCCCTGCCCTTTGTA-3′) and the reverse primer LS266 (5′- GCATTCCCAAACAACTCGACTC-3). The resulting sequences were compared against C. parapsilosis sequences in the NCBI database using BLAST.

Antifungal susceptibility testing

In vitro antifungal susceptibility testing was performed by broth microdilution technique following the European Committee on Antimicrobial Susceptibility Testing (EUCAST) standard methodology for the following antifungals: fluconazole (FLC; Pfizer Inc., Peapack, NJ, USA), itraconazole (ITR, Chemicalpoint, Deisenhofen, Germany), posaconazole (POS; MSD, Rahway, NJ, USA), and voriconazole (VOR; Pfizer Inc., Peapack, NJ, USA). Minimum inhibitory concentration (MICs) were assessed with a nephelometer (Labsystems Nepheloskan Ascent Microplate Reader Type 750) after 24 h of incubation at 35°C. The endpoint of growth was as a 50% inhibition in comparison with the drug-free control. Reference strain C. parapsilosis ATCC 22019 was used as a quality control.

DNA extraction and genome sequencing

For extraction of genomic DNA, single colonies of each putative strain were isolated and grown overnight on YPD plates at 35°C. Fungal material was transferred to a tube containing acid-washed glass beads and lysis buffer (50 mM Tris, 50 mM sodium EDTA, 3% [wt/vol] sodium dodecyl sulfate [SDS] [pH 8]). The samples were homogenized using a vortex adapter, followed by incubation in a thermomixer. Next, cellular debris was pelleted by centrifugation, and the supernatant was transferred to a new tube. RNAse treatment was performed for 1h at 37°C. Subsequently, 3M sodium acetate was added to the tube, samples were incubated at −20°C, then centrifuged and the supernatant transferred to a new tube. DNA was then precipitated using EtOH and the resulting pellet washed twice. DNA extracts that did not meet the requirement of purity were purified using the Monarch genomic DNA purification kit (New England BioLabs Inc., Frankfurt a. M., Germany). Library preparation and 2x250bp paired end sequencing was performed by GeneWiz using the Illumina MiSeq platform.

Genome sequencing, quality control, alignment, and variant identification

Raw sequence reads were trimmed for quality using TrimGalore v0.6.6 in paired end mode using the default error rate threshold of 0.1. Sequence data quality was assessed before and after trimming using FastQC v0.11.9. Mapping to the CDC317 reference genome (FungiDB, release 50) was performed using bwa-mem2 v2.2.1 ⁴⁰. Duplicates were marked using the MarkDuplicates function of Picard v2.24.0. Alignment and genome coverage statistics were assessed using CollectAlignmentMetrics and CollectWgsMetrics from the Picard tools (v2.24.0). Only samples with a mean genome coverage of 10x and >90% of reads mapping to the CDC317 reference genome were included in downstream analyses. Genome coverage and fraction of mapped reads for each sample are listed in Supplemental File 1. Variant calling was performed using Freebayes v1.3.2-dirty. Variants were filtered using vcffilter v0.2 and the following filter string: “QUAL > 30 & DP > 5 & SAF > 0 & SAR > 0 & RPR > 1 & RPL > 1”. Variants were annotated using snpEFF v5.1 ⁴¹ and a custom database made from the CDC317 reference genome files and using the alternative yeast nuclear and yeast mitochondrial codon tables. For defining features within upstream and downstream regions, a cutoff of 500 bp was used.

SNV-based phylogeny and identification of phylogenetic clusters

For the phylogeny and SNV distance-based analyses, only variant positions with non-zero coverage in all samples and SNVs were included. Variant positions with zero coverage, as well as all insertion or deletion variants were removed from the vcf files. The resulting multi-sample VCF was converted to a phylip alignment with a length of 38,567 positions. Pairwise SNV distances between isolates were calculated using SNP Dist ⁴². A maximum-likelihood (ML) phylogeny was inferred by IQ-TREE 2 v2.2.0.3 ⁴³ using the GTR+ASC model to account for ascertainment bias as constant positions were excluded. The ModelFinder module ⁴⁴ of IQ-TREE was used to confirm that GTR was also the best fitting model based on Akaike and Corrected Akaike Information Criterion. The ERG11 amino acid annotations on the phylogeny are reported relative to the C. parapsilosis reference genome CDC317, except for cDNA position 395, corresponding to amino acid position 132. CDC317 is azole resistant and contains a Y132F substitution relative to the EUCAST control strain ATCC 22019. Given this, samples reporting no variant at position 395 of the CDS were re-coded to a Y132F amino acid substitution and samples recording a 395T>A change relative to CDC317 were reported as not having an amino acid substitution at position 132. The annotated phylogenies were created using ggtree ⁴⁵.

Identification of phylogenetic clusters was performed using TreeCluster v1.0.3 ⁴⁶ using a cut threshold of 0.700. The identified clusters were labeled in accordance with the five major clades described in³⁵. As TreeCluster calculated sufficient intraclade genetic distance to further subdivide clades 2 and 4, clusters belonging to these clades were assigned the IDs of 2a, 2b and 4a-g.

Identification of discriminatory loci for sequence typing

Variant profiles of complete gene coding regions were obtained from vcf files and the number of distinct profiles present among the isolates calculated by hierarchical clustering analysis. Protein coding regions with no discriminatory power were discarded and the top 212 CDS were further evaluated for their discriminatory power in 750 bp sliding windows with a step size of 50 bp. Windows including intronic sequence were excluded for downstream selection.

To ensure that potential sequence typing loci were present in the genomes of all samples as a single copy, loci were screened for the presence of the presence of copy number variations (CNVs) using Control-FREEC v11.6 ⁴⁷, DELLY v1.1.6 ⁴⁸, and CNVator v0.4.1 ⁴⁹. Multiple tools utilizing different approaches for detecting copy number changes, including tools utilizing read depth, GC normalization, and split read mapping, were used as many CNV calling software suffer from false positives. DELLY and CNVator were run with default parameters. As FREEC uses GC normalization which varies by organism, 0.25 was used for the minExpectedGC and 0.45 for maxExpectedGC. Loci that fell within CNVs detected by two or more callers were excluded from further consideration. Several samples from the publicly available data displayed prominent chromosome-end bias ⁵⁰. As this overrepresented coverage at the chromosome ends reduces confidence for CNV calling, these samples were excluded from CNV analysis.

In silico strain typing

Locus sequence was extracted in silico from the vcf files using the consensus function of BCFtools v1.9 ⁵¹ and aligned using MUSCLE v3.8 via the R package msa. As BCFtools is not able to convert heterozygous sites to their corresponding IUPAC base, this was performed after alignment using an in-house R script. Concatenated alignments from the four loci were built into a maximum-likelihood phylogeny using RaxML (v2.0.10) using the GTR substitution model and 100 bootstrap replicates. The resulting phylogeny was visualized using ggtree ⁴⁵.

Hunter’s discriminatory power was calculated using the following equation ⁵²:

PCR-amplification and Sanger sequencing of sequence typing markers

Primers for the selected loci were designed using NCBI Primer Blast and are listed in Table 3, along with the predicted fragment size. PCR was performed using the following cycling parameters: 4 min at 95°C, 40 cycles each consisting of a 45 sec denaturation at 95°C followed by 45 sec of annealing at 57°C for CPAR2_808110 and CPAR2_212940 and 59°C for CPAR2_101470 and CPAR2_108720, and an extension at 72°C for 2 min, and a final 10 min at 72°C. Sanger sequencing was performed by LGC (Berlin, Germany). As our Sanger sequence data was not reported with IUPAC codes representing heterozygous sites, these were manually corrected by manual review of the sequencing chromatogram. We observed cases where the CPAR2_101470 sequencing reaction failed to cover the entire amplicon when using the F primer, missing the final variant which we resolved by performing a second Sanger reaction using the R strand. This affects the differentiation of clusters 2b and 2d.