Echinocandin heteroresistance causes prophylaxis failure and facilitates breakthrough Candida parapsilosis infection

Candida parapsilosis is a frequent cause of micafungin breakthrough bloodstream infections

To analyze breakthrough fungal bloodstream infections (BSIs) during allo-HCT, we performed a cohort study that included 952 patients who underwent a total of 983 allo-HCTs from January 2016 to December 2020 at Memorial Sloan Kettering Cancer Center (MSKCC) (Extended Data Fig. 1). In the majority of allo-HCTs (n = 966), patients received antifungal prophylaxis with intravenous micafungin (100-150 mg/day) starting seven to ten days prior to allo-HCT for a median of 17 days, prior to switching to triazole-based prophylaxis per institutional practices ³⁰. We identified a total of 25 micafungin breakthrough BSIs that occurred within 30 days post transplantation (Extended Data Fig. 1), with a clear predominance of Candida parapsilosis as the etiologic agent in these infections (12 out of 25 cases, Fig. 1a).

• Download figure
• Open in new tab

Fig. 1. Micafungin heteroresistance in C. parapsilosis bloodstream isolates.

(a) Micafungin breakthrough bloodstream infections (BSIs) in MSKCC allo-HCT recipients between day 0 and 30 post transplantation (between Jan. 2016 and Dec. 2020). The “others” category includes four episodes of breakthrough BSI by Candida kefyr, Trichosporon asahii, Capronia munkii, and Saccharomyces cerevisiae. (b) Micafungin E-test results of two representative C. parapsilosis isolates, MSK65 (heteroresistant) and MSK804 (susceptible). (c – h) Population-level responses of MSK65 and MSK804 to micafungin. (c, f) Micafungin population analysis profile (PAP) of isolates (c) MSK65 and (f) MSK804. The % survival was calculated compared to growth on drug-free agar. (d, g) CFU dynamics of total and resistant subpopulations of (d) MSK65 and (g) MSK804 in YPD broth containing 2 µg/mL micafungin. Cultures were plated at indicated time points for enumeration of total (black, circles) and resistant (colored, squares) cells. (e, h) Frequency of micafungin-resistant subpopulation of isolates (e) MSK65 and (h) MSK804 during continuous passage. The frequency was calculated after a 24-hour culture in drug-free YPD broth (Pre-treatment), followed by passaging in YPD broth with 2 µg/mL micafungin (In micafungin) for 48 hours, and a second passage in drug-free YPD broth for 72 hours (Subculture). Bar height: mean; Error bar: standard deviation (n = 3). (i) Graphical illustration of the dynamics of the resistant subpopulation in a micafungin heteroresistant isolate in the presence (+ mica) or absence of micafungin (-mica). Blue cells are susceptible, orange cells are resistant (unstable), and grey cells are those killed by micafungin.

Micafungin heteroresistance is prevalent in isolates from C. parapsilosis breakthrough infections

We were able to retrieve 12 C. parapsilosis bloodstream isolates from 10 allo-HCT recipients on micafungin prophylaxis between day 0 and 30 post transplantation (Fig. 1a). In addition, we retrieved three C. parapsilosis isolates from three patients who received micafungin prophylaxis after day 30 of transplantation (Extended Data Fig. 1). Notably, all 15 isolates were classified susceptible to micafungin by broth microdilution (Table 1) and E-test (i.e., Epsilometer test; Fig. 1b) assays. However, neither of these conventional ASTs had sufficient sensitivity to determine if there was a low frequency subpopulation of replicating cells resistant to micafungin, i.e., heteroresistance. Thus, we re-analyzed all isolates by the population analysis profile (PAP) assay, the gold standard for detecting heteroresistance ⁴ (see Methods). In PAP assays, 10⁶ fungal cells were plated on solid agar plates that contained micafungin at indicated concentrations. In this study, we defined an isolate as heteroresistant to micafungin when it contained resistant cells with a frequency between 0.0001% (i.e., limit of detection) and 50% across three consecutive serially diluted micafungin concentrations that included 4 µg/ml (2x breakpoint concentration, consistent with the European Committee on Antimicrobial Susceptibility Testing guidelines ³¹, Supplementary Table 1).

As an example, isolate MSK65 was identified as heteroresistant since it harbored a minor subpopulation of micafungin-resistant cells that survived at concentrations up to 32 µg/ml (Fig. 1c). The resistant cells rapidly replicated and became dominant during micafungin exposure (Fig. 1d), indicating that they were not persister cells which cannot replicate exponentially in the presence of an antimicrobial drug ³². However, the frequency of resistant cells reverted to the baseline level during subsequent growth in the absence of micafungin, indicating that they were not classical resistant cells (Fig. 1e). In contrast, MSK804 was entirely killed by micafungin exposure without the emergence of a resistant subpopulation (Fig. 1f-h). This isolate was classified as micafungin-susceptible. We confirmed that the heteroresistant phenotypes of blood isolates were consistent across three PAP replicates (Extended Data Fig. 2a and 2b).

Using PAP assays, we found that 11 out of the 15 blood isolates were heteroresistant, with the frequency of resistant subpopulations ranging from 0.008% to 0.8% (Table 1). We further determined the MICs for both the susceptible (MIC-S) and resistant (MIC-R) subpopulations in these 11 heteroresistant isolates (see Methods for definition of MIC-S and MIC-R). Due to the predominance of susceptible cells, the clinical MIC values obtained from broth microdilution aligned well with the respective MIC-S values (Table 1). However, the MIC-R values of these isolates ranged from 32 to 64 µg/ml, indicating an increase of 32- to 64-fold when compared to their corresponding MIC-S values. In addition to micafungin, these isolates were generally, but not uniformly, heteroresistant to two other echinocandin drugs, caspofungin and anidulafungin (Extended Data Fig. 2c and 2d). Together, these findings demonstrate the phenomenon of heteroresistance (HR) to echinocandins in C. parapsilosis clinical isolates.

• Download figure
• Open in new tab

Fig. 2. Micafungin heteroresistance and risk of breakthrough C. parapsilosis BSIs in allo-HCT recipients.

(a) A flow chart showing steps involved in determining the number of MSKCC allo-HCT recipients (2016 – 2020) with C. parapsilosis intestinal colonization and breakthrough BSIs. (b) Fisher’s Exact test (two-sided) on patients with micafungin-heteroresistant C. parapsilosis intestinal colonization and breakthrough BSI. (c) Comparative genomic analysis between fecal and blood C. parapsilosis isolates. Longitudinal medical data are shown to the left for five allo-HCT recipients with C. parapsilosis breakthrough BSI, including antifungal drug administration (first row), the time and drug susceptibility of bloodstream isolate(s) (second row), the time of fecal sample collection and the drug susceptibility of corresponding fecal isolate (third row). Pairwise popANI (population average nucleotide identity) between blood (star) and fecal (circle) C. parapsilosis isolates were shown in a heatmap to the right of the medical data. At a cutoff of 99.999%, red colors indicate “identical” isolates, and blue colors indicate “distinct” isolates. The grey dashed line on the diagonal indicate self-comparisons between the same isolates (i.e., popANI = 100%). Patient 1 and Patient 2 in this figure were Patient 2 and Patient 3, respectively, in our previous study ²⁴.

Intestinal colonization by a micafungin heteroresistant strain increases the risk of breakthrough bloodstream infection during prophylaxis

A previous study showed that the intestinal tract is a major reservoir for C. parapsilosis which can translocate across the intestinal mucosa into the bloodstream and cause BSI in allo-HCT recipients ²⁴. To investigate whether micafungin-heteroresistant blood isolates originated from the intestinal tract, we recovered 85 fecal C. parapsilosis isolates from 30 patients (5 of whom developed BSI and were investigated in Fig. 1a). Among these isolates, 43 (50.6%) originated from 14 of the patients and were classified as heteroresistant by PAP (Supplementary Table 2). The remaining 16 of the 30 patients were colonized only by susceptible C. parapsilosis isolates, and none of these patients developed breakthrough BSI. Among the 14 patients who harbored heteroresistant isolates in their intestinal tract, five patients developed a breakthrough C. parapsilosis BSI (Fig. 2a). These data establish an association (P = 0.014; two-sided Fisher’s Exact Test) between intestinal colonization with micafungin-heteroresistant C. parapsilosis isolates and ensuing breakthrough BSI (Fig. 2b).

To confirm the gut-to-blood transmission of micafungin-heteroresistant C. parapsilosis, we sequenced the genomes of fecal and blood C. parapsilosis isolates from the five patients with breakthrough BSIs. In each patient, the bloodstream isolate was identical, defined as >99.999% population average nucleotide identity (popANI) ³³, to an intestinal isolate collected from the feces prior to BSI, with both isolates having the same heteroresistant or susceptible phenotype (Fig. 2c).

Interestingly, two patients (Patient 2 and 3) had a mixture of heteroresistant and susceptible isolates in their feces before BSI (Fig. 2c). Longitudinal analysis revealed that in Patient 2, the heteroresistant isolates replaced the susceptible isolates, leading to an ensuing BSI by a heteroresistant isolate. In Patient 3, both susceptible and heteroresistant isolates translocated to the bloodstream independently, resulting in two separate episodes of BSI (Fig. 2c). These conclusions remained unchanged using an even more stringent popANI cutoff of 99.9999%.

Taken together, our analysis demonstrates that micafungin-heteroresistant C. parapsilosis isolates can colonize the intestine and can subsequently translocate to the bloodstream, thereby presenting a risk factor for the development of breakthrough fungal BSI.

In Patient 2, four of five genomes of susceptible isolates and all genomes of heteroresistant isolates were identical within each category, as defined by the popANI metric. Importantly, the susceptible isolates were distinct from the heteroresistant isolates at the same popANI cutoff (Fig. 2c). Similarly, in Patient 3, two separate clusters of susceptible and heteroresistant isolates were observed. These findings reinforce the notion that the heteroresistance phenotype did not emerge from susceptible isolates during the micafungin prophylaxis period in the patients that we examined. To test this further, we examined 16 patients from whom we isolated at least two identical fecal C. parapsilosis isolates at various time points. Micafungin susceptibility phenotypes remained unchanged among these identical isolates in 15 out of 16 patients (Extended Data Fig. 3). In one patient (patient 47), the phenotype of two identical isolates changed from susceptible to heteroresistance, due to a two-fold change in the MIC of the resistant subpopulation in the later isolate.

Association between micafungin heteroresistance phenotype and genomic variations in clinical isolates from three continents

To expand the scope of the study, we included 118 C. parapsilosis clinical isolates from MSKCC (see Methods) and supplemented the collection with an additional 18 isolates from other hospitals in the United States, 73 isolates from 32 different hospitals from 20 provinces in China, 9 isolates from Germany (Göttingen), and one isolate from France (Nantes), yielding a total of 219 isolates (Supplementary Table 2). Using the PAP method, we found the presence of micafungin-heteroresistant C. parapsilosis isolates in hospitals from all three continents (68/136, 11/73, 7/10 micafungin-heteroresistant isolates from the North America, Asia, and Europe, respectively). This finding highlights the global prevalence of this emerging and clinically significant phenotype.

To understand the genomic underpinnings of micafungin heteroresistance, we sequenced the genomes of these 219 isolates and profiled the genomic variations (GATK, see Methods).

Chromosome-level variations such as aneuploidy have been reported in antifungal heteroresistance in C. parapsilosis and other human fungal pathogens^34–36. We identified 12 out of the 219 isolates (5.5%, eight heteroresistant and four susceptible isolates) with different patterns of trisomy (seven, four, and one isolates have an extra copy of chromosome VI, III, and V, respectively, see Extended Data Fig. 4). Although 8 out of 12 isolates with aneuploidy were micafungin-heteroresistant, the study may lack power to confirm or reject a correlation (two-sided Fisher’s exact test, P = 0.066). An even larger collection of C. parapsilosis isolates would help determine further a relationship between aneuploidy and micafungin heteroresistance. Notably, most heteroresistant isolates (78 of 86) did not exhibit aneuploidy.

We next characterized copy number variations (CNVs), which have been linked with azole tolerance ³⁷, across the 219 C. parapsilosis genomes. We refined a previous algorithm on CNV characterization ³⁸ and defined the copy number of each ORF by two features (see Methods and Extended Data Fig. 5 for illustrative examples): the first feature describes whether an ORF presents a full or a partial amplification or deletion (CNV_categorical or CNV_cat), and the second feature quantifies the copy number of amplified or deleted regions (CNV_quantitative or CNV_quant). With this approach, we identified 3,317 ORFs that contain amplified or deleted regions in at least one isolate in this dataset (Supplementary Table 3). Among the 2,774 ORFs with amplifications, only 53 show overlaps with inverted repeats of at least 100 bp, suggesting that these CNVs were unlikely to be technical artifacts (Supplementary Table 4).

We identified 43,531 single nucleotide variants (SNVs) and 3,722 small insertion/deletions (INDELs) among the 219 isolates (Supplementary Table 5). These variations served as the basis for constructing a phylogenetic tree (see Methods; Fig. 3). Notably, the majority of heteroresistant isolates cluster within one clade, establishing a robust correlation between the heteroresistance phenotype and phylogeny (Moran’s I = 0.30, P < 0.0001).

• Download figure
• Open in new tab

Fig. 3. Phylogenetic tree of 219 C. parapsilosis isolates collected from three continents. The outer circle indicates the micafungin susceptibility phenotype of each isolate.

We next sought to identify variations that are associated with micafungin heteroresistance among the profiled SNVs, INDELs, and CNVs. We first encoded these variations as a mixture of continuous (CNV_quant) and categorical (SNVs, INDELs, CNV_cat) features and then merged redundant features if they correlated perfectly with each other across all 219 genomes (Methods). This resulted in a non-redundant set of 6,806 features (Supplementary Table 6). Using logistic regression (see Methods), we identified 453 features that were significantly associated with the heteroresistance phenotype after multi-test correction (Fig. 4a, Supplementary Table 7). We confirmed that the top five features with both positive and negative associations did not overlap with the inverted repeats identified above (Supplementary Table 4). The most significant feature (Feature 4376) was a homozygous alternative allele at the 5’UTR of CPAR2_206130 (P_adj. = 2.4e-15). The Saccharomyces cerevisiae ortholog of this gene encodes a microtubule-associated protein that plays a role in fungal cell wall organization ³⁹ (Fig. 4b). The second most significant feature (Feature 2154) was the full deletion of CPAR2_600440 (P = 2.9e-15), which encodes a putative glycosylphosphatidylinositol (GPI)-anchored protein involved in cell wall integrity and cell-cell interactions ⁴⁰ (Fig. 4b). Both features remained the top-ranking ones when we examined the associations by LASSO (Least Absolute Shrinkage and Selection Operator; see Methods and Supplementary Table 7).

• Download figure
• Open in new tab

Fig. 4. Top C. parapsilosis genomic features associated with micafungin heteroresistance.

(a) Volcano plot of 6,806 C. parapsilosis genomic features. (b) The top two features positively and negatively associated with micafungin heteroresistance and the annotations of the open reading frames that contain the four features. Distribution of the features were plotted in alignment with the phylogenetic tree. Chromosome locations: Feature 4376: Chr 2; Feature 2154: Chr 6; Feature 3654: Chr 8; Feature 2520: Chr 1. SNV, single nucleotide variation; CNV, copy number variation.

All 453 features with a significant correlation with heteroresistance also correlated with phylogeny (Moran’s I test, P_adj. < 0.05). Among these features, Feature 100 (a homozygous alternative allele within CPAR2_109680, P_adj. = 5.7e-10) was the only one found in an orthologue of the gene encoding the primary drug target. However, this correlation could be due to phylogeny (Extended Data Fig. 6a). We applied a phylogenetic LASSO model to re-evaluate the associations of genomic features with heteroresistance while controlling for their phylogenetic dependences (see Methods). This approach yielded 40 features with non-zero coefficients, where a homozygous alternative variant (Feature 4470) in CPAR2_206800 had the largest positive coefficient (Supplementary Table 8). Although this feature did not completely overlap with the heteroresistant isolates outside the major cluster, it successfully identified the heteroresistance phenotype of three isolates within a phylogenetic region that is dominated by susceptible isolates (Extended Data Fig. 6b). Taken together, the phenotype of micafungin heteroresistance in C. parapsilosis is associated with genomic variations and these associations are strongly, but not exclusively, coupled with phylogeny.

Machine-learning-based diagnosis of micafungin heteroresistance

Given the clinical significance of micafungin heteroresistance, a rapid diagnosis with high accuracy is critical for the rapid administration of effective antifungal drugs. However, current laboratory approaches (i.e., PAP assay) require incubation times up to several days with high material and labor demands. To facilitate the diagnostic process, we developed a machine learning framework (Extended Data Fig. 7) that predicts the micafungin heteroresistance phenotype from the profiled genomic features found in the 219 clinical isolates. Given that multiple isolates from a single patient may share identical genomes (i.e., popANI ≥ 99.999%), separating these isolates into the training and test sets could lead to data leakage and compromise the model evaluation. To address this issue, we removed the duplicates of identical isolates collected from each patient (see Methods, Extended Data Fig. 8), resulting in a total of 160 non-redundant isolates for further analysis. The machine learning framework includes feature standardization, feature selection, hyperparameter optimization, and classification.

In the machine learning framework, we assessed three feature selection methods (logistic regression, LASSO, a neural-network-based approach ENNS ⁴¹) along with three tree-based classifiers: Decision Tree (DT), Random Forest (RF), and Extreme Gradient Boosting (XGBoost). Our approach limited each feature selection method to choose a maximum of ten features to mitigate the risk of overfitting and, more importantly, to reduce the time and effort required to measure the selected features in the laboratory (see Discussion). All model hyperparameters, including the number of features selected, were optimized through a grid search in a 5-fold cross-validation setting. The 160 non-redundant isolates (Supplementary Table 9) were divided into a training set (80%) and a test set (20%), the process of which was repeated 50 times. On average, the combination of logistic regression and XGBoost consistently outperformed the other eight combinations in F1 score (Fig. 5a) and other evaluation metrics (precision, recall, and accuracy, Extended Data Fig. 9a). This combination achieved a median precision of 0.78, a recall of 0.70, an F1 score of 0.74, and an accuracy of 0.84. This strategy can be further simplified by developing SNV-only models, which have equivalent performance to the models developed with SNVs, INDELs and CNVs (Fig. 5a). To confirm that the optimal model captures more than spurious correlations, we randomly permuted the heteroresistance labels to simulate the shuffled data and retrained the model. As expected, the resulting median F1 score never exceeded 0.25 across three repetitions (Fig. 5b). This striking difference strongly indicated that our model effectively captured meaningful genomic features of micafungin heteroresistance.

• Download figure
• Open in new tab

Fig. 5: Evaluation of machine learning model performance in predicting the micafungin heteroresistance phenotype.

F1 score is an evaluation metric that combines precision and recall. Unless specified otherwise, all models selected features from a combined pool of SNV (single nucleotide variant), INDEL (insertion/deletion), and CNV (copy number variation) features. In all boxplots, each dot represents a single train-test split, which was randomly repeated 50 times. (a) Machine learning models constructed by choosing among three feature selection algorithms (LOGIT, Logistic regression; LASSO, Least Absolute Shrinkage and Selection Operator; ENNS, Ensemble Neural Network Selection) and three classifiers (DT, Decision Tree; RF, Random Forest; XGBoost, Extreme Gradient Boosting). (b) Comparison of the best-performing combination (LOGIT+XGBoost) with alterative models. R1, R2, and R3 represent three repeats of retrained LOGIT+XGBoost models using randomly shuffled heteroresistance data. The ANI (average nucleotide identity)-based model quantified genomic similarity between isolates by popANI and assigned a test isolate the same heteroresistance label of its genomically closest isolate in the training set. Phylogeny-based model placed test isolates into a reconstructed phylogeny of training isolates and then performed hidden-state prediction. Two-sided paired-sample t-test: exact P values are provided.

Given the strong association between heteroresistance and phylogeny (Fig. 3), we developed an alternative phylogeny-based predictive approach (see Methods). This approach placed genome sequences of the test isolates into a reference phylogenetic tree created by the training isolates and conducted hidden-state prediction to infer their micafungin heteroresistance labels. Interestingly, the best machine learning model achieved a slightly higher median F1 score (0.74) than the phylogeny-based approach (0.71, Fig. 5b). Although the difference in performance did not attain statistical significance (paired t-test, P = 0.457), both methods substantially outperformed a popANI-based method (Fig. 5b), which simply assigned a test isolate the same heteroresistance label as its most similar (i.e., highest popANI) isolate in the training set. These findings indicate that the features chosen by the machine learning model exploit the phylogenetic signal of heteroresistance and expand on it to achieve predictive power. The most commonly misclassified isolates, including false positive and false negative predictions (Extended Data Fig. 9b), showed micafungin susceptibility phenotypes inconsistent with their phylogenetic neighbors. Collectively, our proof-of-concept machine learning framework demonstrated its effectiveness in predicting the micafungin heteroresistant phenotype based on a small set of genomic features. This innovative approach could facilitate the identification of antimicrobial heteroresistance in clinical settings.

Inclusion criteria for allo-HCT recipients and Candida parapsilosis isolates at MSKCC

Patients that underwent allo-HCT at Memorial Sloan Kettering Cancer Center (MSKCC) between Jan. 1^st, 2016 to Dec. 31^st, 2020 were examined for micafungin breakthrough fungal bloodstream infections. In this study, we defined micafungin breakthrough bloodstream fungal infections as a positive fungal blood culture if a patient received a minimum of one day of micafungin prophylaxis prior to the blood culture positivity, or if micafungin prophylaxis was discontinued within three days before the blood culture positivity. 15 isolates were retrieved from documented breakthrough BSI in these patients: 12 isolates from 10 patients who developed BSI between transplant day 0 to 30 and three isolates from three patients who developed BSI after day 30. The medical record of non-allo-HCT patients with Candida parapsilosis bloodstream infections from the same study period (2016 to 2020) were also reviewed (Supplementary Table 2). A total of 14 bloodstream C. parapsilosis isolates from these non-allo-HCT patients who developed C. parapsilosis BSI based on the criterion of having a positive blood culture were obtained from clinical micro lab at MSKCC.

Fecal biospecimens from allo-HCT recipients were drawn from a prospective fecal collection biobank, which has been described previously ⁴⁷ (MSKCC IRB protocol 16-834). Fecal samples from a subset of 205 patients (collected from −10 to 30 day of transplantation) had been cultured for viable fungal isolates as previously described ²⁴. Among them, a total of 30 patients were confirmed with intestinal colonization by C. parapsilosis from fecal culture. For the five patients with intestinal colonization of C. parapsilosis before BSI, at least one isolate from each fecal sample was randomly picked for further analyses. For the remaining 25 patients with intestinal colonization but no fungemia, one to two C. parapsilosis isolates from each C. parapsilosis culture positive fecal samples were randomly picked. Detailed fecal isolate – patient information data are listed in Supplementary Table 2. Fecal biospecimens from eight healthy donors ⁴⁸ were also examined for the C. parapsilosis culture positivity. Two fecal samples were identified as positive and two C. parapsilosis isolates from each sample were randomly picked. Taken together, there were a total of 89 fecal C. parapsilosis isolates from MSKCC included in this study. Detailed information of patient demographics, antimicrobial exposure, bloodstream infections, and fecal isolates are available in Supplementary Table 2.

C. parapsilosis isolates from non-MSKCC hospitals

Among the 18 C. parapsilosis isolates from United States other than MSKCC, nine isolates were acquired from David Andes’s lab at University of Wisconsin-Madison. The other nine isolates were acquired from the CDC/FDA AR Bank. 73 Asian isolates were acquired from China Hospital Invasive Fungal Surveillance Net (CHIF-NET). These isolates were collected from 32 hospitals from 20 different provinces in China. Ten European isolates include one from Nantes, France and nine from Göttingen, Germany. The detailed information of the isolates are available in Supplementary Table 2.

Classification of micafungin heteroresistance

In this study, an isolate was classified as heteroresistant if the frequency of resistant cells is at least 0.0001%, but remains lower than 50%, across all three consecutive serial-diluted concentrations which include 4 µg/ml (2x breakpoint concentration, EUCAST ³¹). The possible set of three consecutive concentrations could be 1, 2, and 4 µg/ml, or 2, 4, and 8 µg/ml, or 4, 8, and 16 µg/ml. Additionally, we determined the MIC-S (MIC of the susceptible subpopulation) and MIC-R (MIC of the resistant subpopulation) from PAP curves. MIC-S was defined as the minimum concentration at which survival subpopulation drops below 50%. On the other hand, MIC-R was defined as the minimum concentration at which the survival subpopulation falls below the detection limit (0.0001%) at this concentration and all other concentrations beyond it.

Population analysis profile assay (PAP)

YPD agar (Sigma) plates were made with each antifungal drug at the indicated concentrations. The stocks of micafungin (Biovision), caspofungin (Biovision) and anidulafungin (Biovision) were made by dissolving each drug in molecular biology grade water (Corning). Overnight cultures of all the isolates were grown in 2 mL of YPD broth (RPI) in 8 mL round-bottom tubes (Corning) at 35°C with 200 rpm shaking. Ten-fold serial dilutions of the overnight cultures were made in PBS in 96-well plates (Corning). 5 µL of each dilution, from undiluted to 10^-6, were spotted on the agar plates. Plates were incubated for 48 to 72 hours at 35°C before CFUs were enumerated. To calculate the percentage survival, CFUs from the plates containing antifungals were divided by those from the plates without antifungals.

Time-kill assays

2 mL of YPD broth supplemented with 2 µg/mL micafungin for C. parapsilosis or 0.5 µg/mL micafungin for C. auris, were inoculated with 10⁶ CFU from an overnight culture. These YPD cultures with micafungin were incubated at 35°C with shaking at 200 rpm for up to 72 hours. At each indicated time points, 50 µL aliquots were acquired from the cultures and CFUs were counted by plating serial dilutions on plates with and without micafungin at the same concentrations as the respective broth cultures.

Antifungal resistance stability assay

For antifungal treatments, C. parapsilosis cultures were sub-cultured (1:100) from a drug-free overnight culture into YPD with 2 µg/mL micafungin and grown at 35°C with 200 rpm shaking for 48 hours. Micafungin-free subcultures were performed by diluting the micafungin culture 1:1,000 into fresh YPD media and the cultures were grown for 72 hours. Cultures were diluted in PBS and plated on YPD-agar plates with and without micafungin to determine the proportion of resistant C. parapsilosis cells.

E-tests

Micafungin E-test strips were obtained from Biomerieux (Marcy-l’Étoile, France). The E-tests were performed according to the manufacturer’s instructions. Briefly, 5 colonies from an overnight YPD agar plate were dissolved into sterile PBS. The OD was adjusted to 0.5 McFarland standard. The solution was applied to RPMI 1640 plates with MOPS and 2% Glucose (Remel) with a cotton swab. The E-test strip was laid on the plate and the images were captured after 48-hour incubation at 35°C.

Average nucleotide identity (ANI) and identity between isolates

ANI is a quantitative metric for comparing the genomic similarity between two microbial genomes. We calculated the pairwise population ANI (popANI) values for all 219 C. parapsilosis isolates by the inStrain package ³³ (version 1.6.3). A recommended threshold of 99.999% was used to determine whether two isolates are identical or distinct. However, pairwise comparisons may pose challenges in identifying identical isolates when more than two isolates are compared. For example, the identity between isolate A and C cannot be determined when isolate A and B are identical (popANI ≥ 99.999%) and isolate B and C are identical (popANI ≥ 99.999%), but isolate A and C are not identical (popANI < 99.999%).

To address this issue, we developed a computational algorithm to extract the largest clusters in which the popANI value between any two isolates is at least 99.999% and no additional isolate has popANI ≥ 99.999% with all isolates in the cluster. Within each such cluster, all isolates are identical to each other. To identify these clusters, we first conducted hierarchical clustering of the binarized popANI matrix at the 99.999% threshold. Using the “fcluster” function in the Python SciPy package (version 1.10.1), we systematically reduced the maximum allowed inter-cluster distance until the presence of a cluster where all pairwise popANI values among its isolates are at least 99.999%. If multiple such clusters exist, we selected the largest cluster and assigned it a cluster ID. Subsequently, isolates in this cluster were removed from the popANI matrix, and hierarchical clustering was repeated on the remaining isolates. We iterated these steps until no clusters containing more than one isolate could be identified. The remaining isolates were each treated an autonomous cluster. By dividing all isolates into clusters, we define that any two isolates are identical if they belong to the same cluster; otherwise, they are distinct.

Genomic variation calling, filtering, and annotation

Reads were mapped against C. parapsilosis CDC317 reference genome ⁴⁹ using BWA-mem (version 0.7.17-r1188) with parameters -M -R ‘RG’, where ‘RG’ represents read group information generated for each isolate. Duplicate reads were marked using the GATK (version 4.1.9.0) MarkDuplicates tool as per GATK best practices. SNVs and INDELs were identified using the HaplotypeCaller tool from GATK. The resulted GVCF files were combined using the CombineGVCFs tool and jointly genotyped using the GenotypeGVCFs tool from GATK.

The raw SNVs and INDEls were hard-filtered using the following criteria. Filter settings for SNVs include mapping quality (< 40), variant quality score (< 30), quality by depth (< 2), strand odds ratio (> 3), Fisher strand bias (> 60), mapping quality rank-sum test (< −12.5) and read position rank-sum test (< −8). Filter settings for INDELs include quality by depth (< 2), variant quality score (< 30), Fisher strand bias (> 200) and read position rank-sum test (< −20). Very low coverage regions with a read depth < 10 and a genotype quality < 20 were set to missing on a per-sample basis. Monomorphic (no alternative alleles were called for any of the samples) and multiallelic (all sites at which only alternative alleles are called) variants, clusters of three or more SNPs called in a 10-base pair window, and variants in repetitive regions identified by TandemRepeatFinder⁵⁰ (version 4.09.1) were removed. In addition, we removed SNVs that are within 10 bp of INDELs and allowed only one variant within a cluster of INDELs separated by 5 or fewer base pairs to pass. The remaining variants surviving the filtering process were exported in table format using the GATK VariantsToTable tool. SNVs and INDEls were annotated using the SnpEff tool⁵¹ (version 5.0).

Characterization of copy number variation (CNV)

We characterized CNV for each open reading frame (ORF) across eight nuclear chromosomes. The distribution of read depth within each ORF was computed using BEDTools (version 2.31.0)⁵². This distribution was then transformed into a copy number distribution by dividing the read depth by the average read depth across all diploid chromosomes, which was then multiplied by two. Amplified or deleted regions within an ORF can be identified by analyzing the number and position of peaks within the distribution. A normal diploid ORF should display a single peak at the position of two copies, whereas amplifications or deletions result in shifts in peak positions or the emergence of extra peaks. All peaks were identified from smooth Gaussian kernels fitted to the distribution. For each peak, the percentage contribution to the overall distribution was calculated by dividing the area under the peak (boundaries determined as midpoints between this peak and nearest peaks) by the length of the ORF. Small peaks that contributed less than 20% were deemed data noise and filtered out. The remaining peak contributions were then renormalized.

Each CNV was characterized by two variables: a categorical one indicating the type of CNV (normal, full or partial amplification, and full or partial deletion), and a continuous one quantifying the copy number of the amplified or deleted regions. When only one peak remains, the CNV type was categorized as “normal” if the peak centers at two copies, “full amplification” if it centers at a copy number greater than two, or “full deletion” if it centers at a copy number less than two. In cases where two peaks remain and one of them centers at two copies, the CNV type was categorized as “partial amplification” if the other peak centers at a copy number greater than two, or “partial deletion” if it centers at a copy number less than two. If more than two peaks remain and one of them centers at two copies, the amplified or deleted region was determined by the peak with the largest percent contribution among all peaks whose copy number is not two. If two or more peaks remain but none of them center at two copies, both variables were determined based on the location of the largest peak. In the extreme scenario where all peak contributions were below 20% before filtering, these unfiltered small peaks were grouped. The quantitative copy number of this grouped peak was calculated as a weighted average of the copy numbers of these unfiltered peaks.

Identification of C. parapsilosis isolates with aneuploidy

Aneuploidy was identified through the presence of multiple peaks in the distribution of average ORF (open reading frame) coverage across the entire genome. Among the 219 isolates, 15 exhibited two peaks and none exhibited more than two. In isolates with two peaks, the center of the major peak corresponded to the average read depth mapped to normal diploid chromosomes, while the minor peak represented chromosome(s) with potential aneuploidy. The copy number of chromosomes with potential aneuploidy was calculated as the ratio of the minor peak’s read depth to the major peak’s, which was then multiplied by two. We identified 12 isolates with aneuploidy that met the criteria of the copy number ≥ 2.5 or ≤ 1.5.

Identification of inverted repeats

An inverted repeat consists of two arms of similar DNA sequences with one single stranded sequence followed downstream by its reverse complement. The inverted repeats in the C. parapsilosis reference genome CDC317 were identified using Inverted Repeats Finder (version 3.08). The input parameters include Match = 2, Mismatch = 3, Delta = 5, PM = 80, PI = 10, Minscore = 20, MaxLength = 100000, MaxLoop = 1000000.

Feature encoding and merging

SNVs, INDELs, and CNV types are categorical variables and were encoded by one-hot encoding. Specifically, each SNV or INDEL was encoded by two dummy variables: a binary variable indicating whether this variant is a heterozygous allele and a binary variable indicating whether this variant is a homozygous alternative allele. Each CNV type variable was encoded by four dummy variables, each of which indicates whether the type is “full amplification”, “partial amplification”, “full deletion”, or “partial deletion”. No transformation was needed for the quantitative CNV variable.

The total number of variables, including both dummy variables created for SNVs, INDELs, CNV types and unaltered quantitative CNV variables, is 59,881. These variables were termed “raw features”. Additionally, we merged “raw features” that shared identical values across the 219 isolates. These variants that perfectly correlate with each other across multiple genomic loci suggest linkage equilibrium. A total of 6,806 merged features were identified.

Phylogenetic tree construction

The filtered VCF file containing genomic variants served as input for the RRHS tool⁵³ (version, 1.0.0.2). This tool generated haploid FASTA sequences of equal length for each isolate (parameter -n 1 generates one sequence per isolate). Subsequently, the phylogenetic tree was constructed from these sequences using IQ-TREE 2 (version 2.0.7) ⁵⁴. IQ-TREE 2 was executed with the following parameters: -T AUTO (automatically selected threads), -bb 1000 (ultrafast bootstrap approximation ⁵⁵ with 1000 replicates), -m MFP+ASC (use ModelFinder to determine the best fit substitution model with ascertainment bias correction), and -alrt 1000 (1000 bootstrap replicates for SH-like approximate likelihood ration test).

Logistic regression

We used the “Logit()” function implemented in the Python statsmodels package (version 0.14.0) to perform univariate logistic regression, i.e., the association with heteroresistance was assessed for each feature individually. Raw P values were adjusted through the Bonferroni correction. A significance threshold of 0.05 was applied to the adjusted P values.

LASSO and phylogenetic LASSO

The conventional LASSO (Least Absolute Shrinkage and Selection Operator) method determines the coefficient vector β that associates genomic features (matrix X) with heteroresistance label (vector Y) by minimizing ||Y − Xβ||²/2N + α||β||₁, where α is a constant that multiplies the regularization term, and || · ||₁and || · ||₂are L1 and L2 norm respectively. To run LASSO analysis, we called the function “LassoCV” implemented in the Python’s scikit-learn library (version 1.2.2) with CV = 5 (5-fold cross validation) and n_alphas = 1,000 (number of α along the regularization path).

The phylogenetic LASSO model is an extension of conventional LASSO by controlling for variance explained by phylogeny. The technical details of this approach have been described in the Supplementary Information of our previous study for bacterial genomics⁵⁶. Briefly, the matrix X (i.e., genomic features) and vector Y (i.e., heteroresistance label) need to be transformed before applying conventional LASSO: X_new → (T^′T)⁻¹T^′X, Y_new → (T^′T)⁻¹T^′Y, where matrix T is obtained from the phylogenetic variance-covariance matrix G by Cholesky decomposition such as TT^′ = G. The phylogenetic variance-covariance matrix was computed using the “vcv.phylo” function of the ape library⁵⁷ (version 5.7.1).

Machine-learning-based heteroresistance prediction

All machine learning models were developed using data from a training set and subsequently assessed against a test set. As shown above, patient may harbor multiple identical fecal and blood isolates due to repeated sampling and gut-to-blood translocation. If these identical isolates were distributed into both the training and test sets, model performance may be inflated as a result of data leakage. To avoid this potential risk, we systematically excluded duplicate isolates sharing identical genomic sequences (see the section “Average nucleotide identity (ANI) and identity between isolates” for the definition of identical isolates) within each patient. This exclusion process was performed for isolates whose patient information is available.

The Random Forest (RF) and Extreme Gradient Boosting (XGBoost) models were built using a fixed number of 1,000 trees. Other hyperparameters were tuned to maximize the F1 score through grid search in a 5-fold cross-validation framework. The Decision Tree (DT) and RF models shared three parameters, including the number of features chosen by an upstream feature selection algorithm (n_features = 2, 4, 6, 8, 10), the number of features considered when determining the optimal split (max_features = square root (n_features) and n_features), and the maximum depth of each tree (max_depth = 2, 4, 6). The RF models had an extra tunable parameter—the maximum percentage of samples each tree was trained on (max_samples = 0.5, 0.75, 1.0). In the XGBoost models, two unique parameters in addition to n_features and max_depth were optimized: the percentage of features used for building each tree (colsample_bytree = 0.5, 0.75, 1.0) and the learning rate controlling the step size at which model weights were updated (learning_rate = 0.03, 0.1, 0.3).

The implementation of logistic regression and LASSO in Python has been described above. In the case of logistic regression-based feature selection, features were ranked based on their P values (from smallest to largest) and the top “n_features” were selected. When LASSO was employed for feature selection, non-zero features were ranked based on their absolute values (from largest to smallest). In situations where the coefficients of all features are zero, we randomly selected “n_features” features from the entire set of features. When the number of non-zero features was fewer than “n_features”, all non-zero features were selected. The ENNS algorithm has been implemented previously ⁴¹, see code at https://github.com/KaixuYang/ENNS.

The DT and RF models were constructed using the DecisionTreeClassifier and RandomForestClassifier functions from Python’s scikit-learn package (version 1.2.2), respectively. The XGBoost models were developed using the XGBClassifier function from the Python’s XGBoost package (version 1.7.6). We used four metrics, including precision, recall, F1 score (= 2*precision*recall / (precision + recall)), and accuracy, to evaluate the performance of machine learning models.

Phylogeny-based heteroresistance prediction

We used the PICRUSt2⁵⁸ pipeline (version 2.5.2) to predict heteroresistance phenotype. For each train-test split, we began by performing sequence placement using the “place_seqs.py” script. This step involved placing haploid sequences from test isolates into a phylogenetic tree created from training isolates. The custom input files provided to the script “place_seqs.py” included aligned haploid sequences of training isolates (generated by RRHS⁵³, version 1.0.0.2), a phylogenetic tree of the training isolates and a model file specifying optimal tree parameters (generated by RAxML-NG⁵⁹, version 1.2.0), as well as a hidden Markov model of aligned training sequences (generated using hmmbuild from hmmer.org, version 3.1b2). The output of the “place_seqs.py” script is a tree with placed test isolates. With this tree, we next ran hidden-state prediction using the “hsp.py” script to infer their heteroresistance phenotype.

Data visualization

The data obtained from population analysis profile (PAP), time-kill, and antifungal susceptibility assays were graphed by GraphPad prism (San Diego, version 9). All other data were generated and visualized using Python’s matplotlib package (version 3.7.1). All figures underwent further editing in Adobe Illustrator (version 2021). Box plots, swarm plots, and heatmaps were generated using the corresponding functions within Python’s seaborn package (version 0.12.2). In all boxplots, the central line represents the median, the box limits denote the first and third quartiles respectively, and the whiskers extend to the smallest and largest values or at most to 1.5× the interquartile range, whichever is smaller.

Other computational tools

In general, Python version 3.8.5 and R version 4.3.1 were used in this study for computational analyses. Moran’s I test for phylogenetic correlation was performed using “Moran.I” function from the “ape” package (version 5.7.1). Phylogenetic trees were plotted and annotated using the package “ggtree” (version 3.8.0). Fisher’s exact test and paired-sample t-test were performed using “fisher_exact” and “ttest_rel” from the SciPy (version 1.10.1) stats module respectively.