1. Computational and Systems Biology
  2. Microbiology and Infectious Disease

Candida albicans exhibits heterogeneous and adaptive cytoprotective responses to antifungal compounds

  1. Vanessa Dumeaux
  2. Samira Massahi
  3. Van Bettauer
  4. Austin Mottola
  5. Anna Dukovny
  6. Sanny Singh Khurdia
  7. Anna Carolina Borges Pereira Costa
  8. Raha Parvizi Omran
  9. Shawn Simpson
  10. Jinglin Lucy Xie
  11. Malcolm Whiteway
  12. Judith Berman
  13. Michael T Hallett  Is a corresponding author
  1. Department of Anatomy and Cell Biology, Western University, Canada
  2. Department of Biology, Concordia University, Canada
  3. Department of Computer Science and Software Engineering, Concordia University, Canada
  4. Shmunis School of Biomedical and Cancer Research, The George S. Wise Faculty of Life Sciences, Tel Aviv University, Israel
  5. Department of Chemical and Systems Biology, Stanford University, United States
  6. Department of Biochemistry, Western University, Canada
https://doi.org/10.7554/eLife.81406
Share this article
Cite this article
  1. Vanessa Dumeaux
  2. Samira Massahi
  3. Van Bettauer
  4. Austin Mottola
  5. Anna Dukovny
  6. Sanny Singh Khurdia
  7. Anna Carolina Borges Pereira Costa
  8. Raha Parvizi Omran
  9. Shawn Simpson
  10. Jinglin Lucy Xie
  11. Malcolm Whiteway
  12. Judith Berman
  13. Michael T Hallett
(2023)
Candida albicans exhibits heterogeneous and adaptive cytoprotective responses to antifungal compounds
eLife 12:e81406.
https://doi.org/10.7554/eLife.81406
  2. Download BibTeX
  3. Download .RIS

Abstract

Candida albicans, an opportunistic human pathogen, poses a significant threat to human health and is associated with significant socio-economic burden. Current antifungal treatments fail, at least in part, because C. albicans can initiate a strong drug tolerance response that allows some cells to grow at drug concentrations above their minimal inhibitory concentration. To better characterize this cytoprotective tolerance program at the molecular single-cell level, we used a nanoliter droplet-based transcriptomics platform to profile thousands of individual fungal cells and establish their subpopulation characteristics in the absence and presence of antifungal drugs. Profiles of untreated cells exhibit heterogeneous expression that correlates with cell cycle stage with distinct metabolic and stress responses. At 2 days post-fluconazole exposure (a time when tolerance is measurable), surviving cells bifurcate into two major subpopulations: one characterized by the upregulation of genes encoding ribosomal proteins, rRNA processing machinery, and mitochondrial cellular respiration capacity, termed the Ribo-dominant (Rd) state; and the other enriched for genes encoding stress responses and related processes, termed the Stress-dominant (Sd) state. This bifurcation persists at 3 and 6 days post-treatment. We provide evidence that the ribosome assembly stress response (RASTR) is activated in these subpopulations and may facilitate cell survival.

Editor's evaluation

The valuable study by Dumeaux et al. examines the transcriptional response to antifungal treatment in the major opportunistic human fungal pathogen Candida albicans. Using solid methodology, including a novel droplet-based single cell transcriptomics platform, the authors report that fungal cells exhibit heterogeneity in their transcriptional response to antifungal drug treatment. The ability to study the trajectories of individual cells in a high-throughput manner provides a novel perspective on studying the emergence of drug tolerance and resistance in fungal pathogens.

https://doi.org/10.7554/eLife.81406.sa0

eLife digest

Many drugs currently used to treat fungal diseases are becoming less effective. This is partly due to the rise of antifungal resistance, where certain fungal cells acquire mutations that enable them to thrive and proliferate despite the medication. Antifungal tolerance also contributes to this problem, wherein certain cells can continue to grow and multiply, while other – genetically identical ones – cannot. This variability is partly due to differences in gene expression within the cells. The specific nature of these differences has remained elusive, mainly because their study requires the use of expensive and challenging single-cell technologies.

To address this challenge, Dumeaux et al. adapted an existing technique to perform single-cell transcriptomics in the pathogenic yeast Candida albicans. Their approach was cost effective and made it possible to examine the gene expression in thousands of individual cells within a population that had either been treated with antifungal drugs or were left untreated.

After two to three days following exposure to the antifungal treatment, C. albicans cells commonly exhibited one of two states: one subgroup, the ‘Ribo-dominant’ cells, predominantly expressed genes for ribosomal proteins, while the other group, the ‘Stress-dominant’ cells, upregulated their expression of stress-response genes. This suggests that drug tolerance may be related to different gene expression patterns in growing cell subpopulations compared with non-growing subpopulations.

The findings also indicate that the so-called ‘ribosome assembly stress response’ known to help baker’s yeast cells to survive, might also aid C. albicans in surviving exposure to antifungal treatments.

The innovative use of single-cell transcriptomics in this study could be applied to other species of fungi to study differences in cell communication under diverse growth conditions. Moreover, the unique gene expression patterns in C. albicans identified by Dumeaux et al. may help to design new antifungal treatments that target pathways linked to drug resistance.

Introduction

Candida albicans is one of the most prevalent human fungal pathogens (Benedict et al., 2019; Kullberg and Arendrup, 2015; Pappas et al., 2018). Systemic C. albicans infections are the second most common cause of mortality from infectious diseases in extremely premature infants (15–20% mortality), and the fourth most common cause of nosocomial bloodstream infections (30–50% mortality) (Benjamin et al., 2010; Brown and Netea, 2012; Pfaller et al., 2019). The frequency of drug resistance is far lower than the frequency of treatment failures, such that most infections that fail to respond to antifungal drug treatments are susceptible to the drug. This clinical persistence can result from heterogeneity in antifungal drug responses, variations in host immune status, and the inability of drugs to reach their fungal target sites (Delarze and Sanglard, 2015; Wuyts et al., 2018).

A resistant fungal isolate can thrive in the presence of an antifungal drug at concentrations exceeding the minimum inhibitory concentration (MIC) of that species, with detectable growth within 24 hr of drug exposure. For a tolerant fungal isolate, a proportion of cells (typically 5–95% of the population) can grow, albeit slowly, at concentrations surpassing the population average MIC (Berman and Krysan, 2020); tolerance is detectable only 48 hr after drug exposure, which is longer than most clinical assays are performed (Espinel-Ingroff et al., 2017). Approaches for predicting patient response to an infection, based solely on the average response of cells to a drug challenge across a population after 24 hr, can overlook such subpopulations with behaviors that have the potential to significantly impact the evolution of drug resistance (Altschuler and Wu, 2010).

Resistance usually results from mutations that directly influence interaction of the drug with its target (Berman and Krysan, 2020), while tolerance is thought to stem from phenotypic heterogeneity or cell-to-cell variations in phenotypic responses within an isogenic cell population. Microbial phenotypic heterogeneity can arise from mechanisms including stochastic or periodic gene expression, protein stability, cell age, cell-cell interactions, chromatin modifications, and genomic neoplasticity (Ackermann, 2015). Heterogeneity in microbial populations can confer benefits like bet-hedging, rapid metabolic shifts, division of labor, and resource sharing (Ackermann, 2015). It is important to note the distinction between antibacterial and antifungal drug tolerance. For bacteria, tolerance is defined by the duration of cell survival following periodic exposure to a bactericidal drug (Brauner et al., 2016). By contrast, for fungi, tolerance is characterized by the proportion of growth in supra-MIC drug concentrations relative to growth without the drug and is best characterized in C. albicans treated with fluconazole (FCZ). The underlying mechanisms and implications of antifungal drug tolerance differ from those of bacterial heteroresistance, in which a small subpopulation exhibits growth in bacteriocidal drugs (Andersson et al., 2019). In this study, the focus is on antifungal tolerance to FCZ, a well-known fungistatic drug.

Previous studies have characterized key aspects of drug tolerance in C. albicans. Tolerance is largely drug concentration-independent and high tolerance correlates with poor clinical outcomes (Astvad et al., 2018; Levinson et al., 2021; Rosenberg et al., 2018). Tolerance increases the effective population size (Cowen and Lindquist, 2005; Vincent et al., 2016). The emergence of higher tolerance cannot be attributed solely to the accumulation of adaptive point mutations (Wertheimer et al., 2016), as selected tolerant cells and their original parent give rise to mixed populations of tolerant and non-tolerant cells in similar proportions (Rosenberg et al., 2018). Adjuvant drugs used alongside the common fungistatic FCZ not only inhibit tolerance but prevent the evolution of resistance because these drug combinations are fungicidal. Such adjuvant drugs can modulate the tolerance response by targeting Hsp90, calcineurin, TOR, PKC, and sphingolipid biosynthesis (Cowen and Lindquist, 2005; Karababa et al., 2006; Rosenberg et al., 2018; Sanglard et al., 2003; Vincent et al., 2016). The molecular mechanisms coordinating tolerance regulation across cellular pathways are not yet fully understood, although stress responses appear to be central (Cowen and Lindquist, 2005; Rosenberg et al., 2018). Several cellular processes are essential for tolerance, including ergosterol biosynthesis and proteasome function (reviewed in Berman and Krysan, 2020).

Most previous studies of drug tolerance, ranging from classic growth assays to omics profiling, have treated surviving fungal subpopulations as a homogenous whole, mainly due to a lack of technology to efficiently study large numbers of individual fungal cells. While single-cell (sc) transcriptomic assays have been applied extensively to mammalian systems, their use in fungal contexts remains limited due to the challenges of disrupting the rigid cell wall, lysing the membrane, and the small overall amount of RNA per cell. In the model yeast Saccharomyces cerevisiae, several limited microfluidic or barcoding sc studies examined approximately a hundred cells (Gasch et al., 2017; Nadal-Ribelles et al., 2019b; Nadal-Ribelles et al., 2019a; Urbonaite et al., 2021), and two high-throughput fungal transcriptome studies profiled ~40K cells (Jackson et al., 2020; Jariani et al., 2020). In our preliminary pre-printed study (Bettauer et al., 2020), we profiled C. albicans using a fungal nanoliter droplet-based assay (DROP-seq), modified from the original system presented by Macosko et al., 2015. This approach overcomes technical challenges that arise in fungal settings, providing a flexible, cost-effective solution. Here, we extend the data and analysis from Bettauer et al. to better understand subpopulation-specific responses to drug stress. Since the original sc study (Bettauer et al., 2020), a second C. albicans sc study was performed to investigate the initial response (within 3 hr) to high FCZ concentrations (Dohn et al., 2022). Dohn and colleagues observed increased expression of ergosterol genes, an increase in cell cycle arrest genes, and a transient increase in stress response gene expression at 1.5 hr. The work also identifies interesting variability in the temporal response to acute drug exposure but does not address molecular heterogeneity in the subsequent emergence from cell cycle arrest that can only be observed as tolerance after longer times (days 2–3), the focus of this study.

Here, we use the fungal DROP-seq system (Bettauer et al., 2020) to explore phenotypic heterogeneity in the C. albicans response to antifungal drugs. We profiled the transcriptomes of thousands of individual cells from C. albicans populations that were either untreated (UT) or exposed to antifungal compounds with a focus on FCZ across several days. This data is integrated with bulk DNA sequencing and fluorescence microscopy to provide in-depth analyses that focus on subpopulation composition and phenotypic heterogeneity across isogenic cell populations. We detect a heterogeneous response to FCZ within isogenic populations, supporting the idea that different cells exhibit distinct survival strategies, including some with increased expression of genes related to drug tolerance. This study underscores molecular events that may lead to drug tolerance and that hold potential for future therapeutic targeting.

Results

C. albicans exhibits drug tolerance 2 days after exposure to FCZ

To uncover molecular events associated with the emergence of drug tolerance in C. albicans, we performed a series of disk diffusion assays, focusing on a 6-day time series with FCZ, detailed in Appendix 1. When assayed at 48 hr, both disk assays and broth microdilution assays report on tolerance (Rosenberg et al., 2018). Drug tolerance increased significantly between days 1 and 2 (p<0.001, Kruskal-Wallis χ2 test) and increased slightly at day 3, after which the tolerance level did not change. Thus, tolerant cells are present in all cultures starting at day 2. Since these populations originate from a single isogenic colony, we assume that tolerant and non-tolerant cells in the same culture differentially express pathways and processes relevant to their ability to grow (or not grow) in the presence of the drug. We do not observe widespread drug resistance at any time point.

An optimized sc profiling assay to explore drug tolerance in C. albicans

Although sc profiling with a commercial system is feasible in S. cerevisiae (Jackson et al., 2020), specific aspects of fungal biology motivated us to develop a low-cost alternative tailored for fungi, specifically C. albicans. We first optimized the removal of the cell wall, as well as the time and concentration parameters to fix the transcriptome (Materials and methods, ‘Strains, media, and drug treatment’, ‘Spheroplasts’). We then built a nanoliter droplet-based system modified from Macosko et al., 2015, using in-house components as described previously (Booeshaghi et al., 2019), to reduce the device and assay costs. The profiles reported here combine our preliminary effort (Bettauer et al., 2020) with additional data and analyses to provide increased power to examine the technical and biological efficacy of our system.

C. albicans populations were grown in rich media alone (UT) or with an antifungal compound: FCZ (1 µg/ml), caspofungin (CSP; 1 ng/ml), or rapamycin (RAPA; 0.5 ng/ml) (Materials and methods, ‘Strains, media, and drug treatment’). UT samples were collected during the logarithmic growth phase and treated samples were collected at days 2 and 3 post-drug exposure. This is the period when FCZ tolerance becomes evident. To explore whether subpopulations observed at days 2 and 3 persisted at later time points, we resuspended the FCZ day 3 population in fresh medium and profiled the samples on day 6 after 3 days of growth without drug (Figure 1A). The vast majority of cells were in the yeast white morphology with less than 0.2% of cells exhibiting a filamentous morphology (hyphae or pseudohyphae).

Figure 1 with 2 supplements see all
Download asset Open asset
Experimental design and initial single-cell profiling.

(A) The time series experiment begins with three replicates of untreated (UT) cells followed by profiling of rapamycin (RAPA), caspofungin (CSP), and fluconazole (FCZ, 2 replicates at day 2, and 3 replicates at day 3). After 3 days in FCZ, cells were transferred to YPD; recovered cells were profiled at day 6 (i.e., 3 days after resuspension). (B) Bar plot (left) depicts the number of high-quality cells per sample. Violin plots (right) the distribution in the number of reads assigned to each cell.

After processing with the DROP-seq device, samples were sequenced following the original protocol (Macosko et al., 2015) but with cell concentration and PCR cycle numbers optimized for C. albicans (Materials and methods, ‘Cell preparation for sc profiling’). Sc sequence data was processed using a reference index that covers the spliced transcriptome and the Alevin-Fry package (He et al., 2021) (Materials and methods, ‘Quality control, basic processing, and normalization of the sc profiles’). Gene and cell quality control are challenging exercises in all sc profiling efforts (Svensson, 2019), and especially in fungi because of the small amount of RNA per cell, especially under the stress of growth in antifungal drugs (Jariani et al., 2020). A series of quality control procedures were used to estimate gene/cell expression levels (Figure 1—figure supplement 1A, Materials and methods, ‘Quality control, basic processing, and normalization of the sc profiles’).The pipeline identified 18,854 high-quality cells across the 11 drug/time point conditions with an average of 1714 cells per sample. On average, 184 transcripts were identified in each cell, however there is large dispersion in the right tail representing many cells with significantly more transcripts (max. 1984) (Figure 1B). On average, these transcripts arose from 94 unique genes per cell, again with large right tail dispersion (max. 825) (Figure 1—figure supplement 1B and C). Since theoretical results highlight the importance of many cells over the number of identified genes per cell (Zhang et al., 2020), we reasoned that inclusion of the sparse cells (left tail) would strengthen the analyses and help identify large subpopulations with strong differentially expressed transcriptional programs across the different treatments. Moreover, although the gene by cell count matrix was sparse, there was very high concordance between FCZ pseudo-bulk profiles (i.e., aggregated sc counts) at day 2 and day 3 (R=0.82; Figure 1—figure supplement 2A), indicating that the assay robustly quantifies the expression of genes across different batches.

To further investigate the robustness of the assay, we performed bulk RNA-sequencing of FCZ-treated cells at day 2 post-treatment (Materials and methods, ‘Bulk transcriptomics’) and compared this bulk profile with the pseudo-bulk derived from mapping sc reads to the reference genome but ignoring barcodes (unfiltered ‘pseudo-bulk’ profiles; Materials and methods, ‘Construction of pseudo-bulk profiles’). The comparison identified 6071 genes with only 172 genes not detected in one or more of the pseudo-bulk profiles. We note that missing genes were mostly lowly expressed genes in the bulk profile (Figure 1—figure supplement 2B). Moreover, day 2 and day 3 pseudo-bulk replicates were significantly correlated with the bulk RNA-sequencing (Figure 1—figure supplement 2C). This strongly suggests that the DROP-seq-derived profiles sample the C. albicans transcriptome, capture true biological signals, and primarily miss transcripts expressed at lower levels.

In isogenic UT cells, differences in cell transcriptional profiles highlight metabolic and stress responses coupled with cell cycle checkpoints

As described above, mid-log phase cells grown under standard conditions were collected for sc profiling. In addition to sc transcriptomes, ‘bulk’ DNA-sequencing profiles were generated to verify strain isogenicity (Materials and methods, ‘Whole-genome DNA-sequencing’).

To identify the main sources of cell-to-cell variability in UT cells (N=5062 cells), we compiled gene signatures for biological processes and responses likely to play a role in microbial phenotypic heterogeneity based on previous transcriptional profiling studies. Gene signatures often consist of genes that are differentially expressed when a specific process is activated (or repressed) compared to wildtype cells. Here, the signatures included processes such as cell cycle, stress responses (both specific and general), the TCA cycle, and metabolic pathways (such as glycolysis) (Figure 2—source data 1, Materials and methods, ‘Cell clustering, trajectory, and signature analyses’). Most of these expression signatures were derived in the context of bulk transcriptional studies either directly from C. albicans or from other fungi. Where necessary, we identified orthologs of relevant genes in C. albicans (Balakrishnan et al., 2012), while considering regulatory differences between species (Johnson, 2017; Lavoie et al., 2010).

The score of each signature was measured in each UT cell profile (Materials and methods, ‘Cell clustering, trajectory, and signature analyses’); the most variable signatures are displayed in Figure 2A and B. One large group of cells exhibited elevated expression of genes involved in both M phase and the heat-shock response (cell indices >3000), consistent with a study linking these two processes (Senn et al., 2012). We also found that a small group of cells with very high expression of heat-shock response exhibited the lowest expression of M and S phase, glycolysis, and RP-coding-related genes (cell indices ~2700–3000). These few cells likely experienced cell cycle arrest due to high levels of cellular stress. Most other M phase cells with relatively high expression of heat-shock protein-coding genes (indices 3000–5000) also exhibited high expression of glycolysis-related genes. In contrast, cells with the least evidence of M phase expression, sometimes exhibited high S phase gene expression, low expression of the heat-shock signature, and relatively high expression of the oxidative stress signature (cell indices 0–2700, Figure 2A and B). Overall, these results indicate that expression heterogeneity in UT cells is primarily related to different cell cycle phases, distinct stress responses, and metabolic states. Of note, the UT cells were grown in rich media without exposure to any known stressors.

Figure 2 with 1 supplement see all
Download asset Open asset
Cell-to-cell heterogeneity in untreated (UT) cell populations.

(A, B) Expression levels (VISION scores) of curated signatures for individual UT cells. Cell order is the same in all graphs. (C) Scatterplot of cells (dots) based on expression level after imputation. Colors indicate: red, TTR1 expression >1.2 and HSP70 expression <1.2; green, TTR1 expression <1.2 and HSP70 expression >1.2; black, both TTR1 and HSP70 expression >1.2; and gray, expression of both genes was <1.2. Distributions of expression are illustrated in histograms above and to the right and the number of cells in each group is provided in the top right of the figure. (D) Representative fluorescence micrographic image of RFP-tagged TTR1(red) and GFP-tagged HSP70 (green) in a population of isogenic cells showing the mutually exclusive nature of their expression in individual cells. (E) Plot of the mean intensities captured in microscopy images of RFP-tagged TTR1 (red) and GFP-tagged HSP70 (green).

Figure 2—source data 1

Gene signatures related to microbial phenotypic diversity and drug tolerance curated from the literature.

The most variable signatures identified in each analysis are highlighted – dark colors: final selection; light colors: signatures passed the threshold but not selected after filtration (redundancy; non-concordant gene expression; number of most variable genes ≤1).

https://cdn.elifesciences.org/articles/81406/elife-81406-fig2-data1-v2.xlsx
Figure 2—source data 2

GFP and RFP oligo primers.

https://cdn.elifesciences.org/articles/81406/elife-81406-fig2-data2-v2.xlsx

To further investigate the hypothesis that UT cells differentially express genes involved in specific stress responses, we chose gene pairs predicted to have mutually exclusive expression in any given cell based on the sc transcriptomics profiles; for example, heat-shock protein 70 (HSP70) and dithiol glutaredoxin (TTR1), which is involved in the oxidative stress response (Figure 2C; McNemar test, p-value <0.001). We constructed a dual fluorescent reporter strain expressing GFP-labeled Hsp70 and RFP-labeled Ttr1 (Materials and methods, ‘Cell imaging’). During growth in rich medium, fluorescence microscopy revealed a notable level of expression with only one of the two markers detectable in a given cell as predicted from the sc profiles (Figure 2D, E, Figure 2—figure supplement 1, McNemar test, p-value <0.001). Thus, the distinct transcriptional stress responses observed in the UT population data are not due to stress responses during sc profiling. Instead, they most likely represent the cell cycle stage and/or the metabolic condition of the cell, emphasizing the significance of metabolic- and stress-sensitive phases in the progression of the cell cycle. These results align with several previous reports that associate cell cycle phase with the expression of stress response and metabolism-related genes (Brauer et al., 2008; Chiu et al., 2011; Hossain et al., 2021; Senn et al., 2012).

Cells display distinct survival responses to FCZ

We next investigated the C. albicans response to antifungal compounds using fungistatic FCZ (1 µg/ml, 1–2 × MIC50, Figure 3—figure supplement 1A), and fungicidal CSP (1 ng/ml, <0.03 × MIC). These concentrations were chosen to ensure a sufficient number of survivors. FCZ targets Erg11p, which encodes an enzyme central to ergosterol biosynthesis and membrane integrity (Odds et al., 2003; Thamban Chandrika et al., 2018; Wertheimer et al., 2016). At these concentrations, FCZ treatment slows growth relative to UT controls in the first days after exposure (Figure 3—figure supplement 1B). Since tolerance to FCZ should be evident after 2 days of incubation in the presence of the drug (Gerstein and Berman, 2020), we profiled our C. albicans cell subpopulations at 2 and 3 days post-exposure using drug disk assays (Appendix 1). Finally, we also investigated the C. albicans response to low doses of RAPA (MIC80 <1 µg/ml, 0.5 ng/ml; 0.0005× MIC) (Cruz et al., 2001). RAPA is known to inhibit tolerance to FCZ (Rosenberg et al., 2018).

We focused on 11,309 good-quality cells captured at 2 or 3 days post-exposure to FCZ, CSP, and RAPA, in addition to the 5062 UT cells. Unsupervised analysis was used to identify similar clusters of cells. The analysis first uses a deep generative neural network (scVI) to transform the high-dimensional expression profiles down to lower dimensions, followed by Leiden clustering to identify groups of cells that have similar expression profiles (Materials and methods, ‘Quality control, basic processing, and normalization of the sc profiles’, ‘Cell clustering, trajectory, and signature analyses’). This revealed five main clusters containing 92% of all N=16,371 cells, in addition to ~20 small clusters. We used Uniform Manifold Approximation and Projection (UMAP) to visualize these patterns in 2D; cells are colored by their cluster assignment (Figure 3A). To investigate the resilience of the cell clusters, we randomly repeated the clustering process 100 times, each time using a random subset consisting of 95% of the cells. Cells originally assigned to four of the main clusters almost always remained in that cluster (Figure 3—figure supplement 1C). A significant portion of cells allocated to the (smaller) cluster 5 (5-purple, N=501 cells, 3%) were categorized together with cluster 3 cells, which is indicative of only subtle gene expression differences and suggests the two clusters are collapsed (Figure 3—figure supplement 1C). The most variable signatures across clusters are displayed in Figure 3C, where color corresponds to the average score (z-score, color bar) across all cells within each cluster. UT cells are primarily found in cluster 3-green and to a lesser extent in cluster 2-pink (Figure 3B). CSP cells were mostly found in cluster 1-darkpink while RAPA cells were mostly found in cluster 2-lightpink (Figure 3B). Biological inferences associated with each of these clusters are further described in Appendix 2. The remaining 3% of cells (N=573) were outliers scattered across ~20 distinct clusters (‘comet’-like cluster-darkblue; Figure 3A). The ‘comet-like’ clusters appear in random directions from the five main clusters. This pattern suggests that the small set of cells in each comet have strong transcriptional similarity, but each such comet is transcriptionally distinct from other comets. The comets were primarily observed in FCZ survivors at 2 days (Figure 3B); an investigation of the relevant biology underlying comets is discussed in Appendix 2.

Figure 3 with 2 supplements see all
Download asset Open asset
Cell profiles after challenges with different antifungal compounds.

(A) Uniform Manifold Approximation and Projection (UMAP) embedding of all cells including untreated (UT) cells and cells treated with fluconazole (FCZ) (days 2 and 3), rapamycin (RAPA) (day 2), and caspofungin (CSP) (day 2). Leiden clustering identified five main clusters and ~20 sparsely populated ‘comets’. (B) Dotplot describing relative size of cluster populations. Rows correspond to clusters and columns correspond to drug_day conditions. Dot diameter is proportional to the fraction of cells from each condition for a given cluster. The dot color is proportional to the fraction of cells from each cluster for a given condition as shown in the color bar. Numbers in parentheses indicate the total cell count in clusters and drug_day conditions. (C) A heatmap depicting the level of activation (VISION z-scores) of different signatures. (D–G) The UMAP embedding from panel (A) but where color depicts: (D) density of FCZ at day 2 cells, (E) density of FCZ at day 3 cells. (F,G) Signature scores for expression of (F) ribosomal protein (RP), (G) and heat-shock stress (HSP) signatures.

Figure 3—source data 1

Differentially expressed genes in fluconazole (FCZ)-treated cells classified in the Ribo-dominant (Rd) cluster 1 (pseudo-bulk samples) compared to FCZ-treated cells in the Stress-dominant (Sd) cluster 4 (pseudo-bulk samples).

baseMean = the average of the normalized counts taken over all pseudo-bulk samples; log2FoldChange = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange estimate; stat = Wald statistic; pvalue = Wald test p-value; padj = Benjamini-Hochberg adjusted p-value.

https://cdn.elifesciences.org/articles/81406/elife-81406-fig3-data1-v2.xlsx
Figure 3—source data 2

Gene ontology terms enriched in genes differentially expressed between fluconazole (FCZ)-treated cells classified in the Ribo-dominant (Rd) cluster 1 (pseudo-bulk samples) compared to FCZ-treated cells in the Stress-dominant (Sd) cluster 4 (pseudo-bulk samples) listed in Figure 3—figure supplement 2B.

https://cdn.elifesciences.org/articles/81406/elife-81406-fig3-data2-v2.xlsx

Interestingly, a clear bifurcation in the FCZ survivor cell population was evident in the sc transcriptional profile by day 2 (Figure 3D), corresponding to clusters 1-darkpink and 4-turquoise of Figure 3A. We term these two distinct states the Ribo-dominant (Rd) state and the Stress-dominant (Sd) state, respectively. The bifurcation between the Rd and Sd states becomes more pronounced by day 3 (Figure 3E), with nearly every surviving FCZ-treated cell appearing in either the Rd or Sd state.

The Rd state is characterized by high RP gene expression, moderate to high expression of GCN4-mediated response genes that activate amino acid biosynthesis, along with low expression of glycolytic and carbohydrate reserve metabolic pathway genes, and an absence of heat-shock or hyperosmotic stress response genes (Figure 3C-G; Figure 3—figure supplement 2A). By contrast, the Sd state is characterized by high expression of heat-shock stress response genes and low expression of RP genes (Figure 3C-G; Figure 3—figure supplement 2A).

In total, 797 genes are differentially expressed between the Rd and the Sd state (pseudo-bulk from cluster FCZ 1-dark pink versus cluster 4-turquoise; DESeq2, FDR <0.1; Materials and methods, ‘DGE analysis’; Figure 3—source data 1 ). Among these genes, 230 are overexpressed in Rd cells, with half of them involved in protein translation (Figure 3—source data 1 and 2; Figure 3—figure supplement 2B). Many highly expressed genes within the Rd response were also involved in rRNA processing and mitochondrial cellular respiration.

By contrast, 567 genes are more highly expressed in the Sd state; this set of genes is enriched for involvement in cell wall organization, cell adhesion, morphology, virulence, filamentous growth, and biofilm formation (Figure 3—source data 1 and 2; Figure 3—figure supplement 2B). Interestingly, Sd cells also highly express genes involved in the unfolded protein response (UPR), such as HSP70 and YHB1, as well as genes that promote drug tolerance, such as HSP90, GZF3, CCH1, HSP21, HSP70, and RIM101 (Delarze et al., 2020; Garnaud et al., 2018; Liu et al., 2015; Mayer et al., 2013; Nagao et al., 2012; Rosenberg et al., 2018).

The sc profiles revealed that isogenic cells that survived FCZ treatment are found in one of two distinct states at days 2 and 3. This highlights the heterogeneity of cellular responses to FCZ treatment, which may be associated with drug tolerance and clarifies that one state (Rd) is enriched in ribosome and translation-related functions while the other state (Sd) is enriched in genes related to genes identified in the stress responses including the UPR and drug tolerance.

Drug tolerance in C. albicans may involve the ribosome assembly stress response to facilitate a transition from the Rd to Sd state

The Rd state, evident on days 2 and 3, is marked by increased expression of RP and rRNA processing genes. Proper ribosomal biogenesis requires a balance between the synthesis of RP and rRNA. An imbalance between RP and rRNA processing can lead to proteotoxic stress and RP aggregation (Tye et al., 2019). This stress and aggregation activate the ribosome assembly stress response (RASTR) in yeast, which involves heat-shock transcription factor Hsf1 (Albert et al., 2019). Hsf1 subsequently upregulates HSP90 and other protein folding-related genes, including HSP70, which are upregulated in the Sd state cells. In this context, we propose that the drug tolerance response in C. albicans may involve RP aggregation (in Rd cells), which activates RASTR, and induces the expression of Hsf1 (in Sd cells). This process might allow cells to transition from the Rd to Sd state, where genes implicated in drug tolerance are expressed.

To test this hypothesis, we collected a list of C. albicans orthologs of S. cerevisiae RASTR genes (gene list supplied by B Albert, Table 1A) together with a list of C. albicans genes which are targets of Hsf1 (Leach et al., 2016; Table 1B). In S. cerevisiae, the characteristic RASTR gene expression profile includes decreased expression of several RP genes and increased expression of Hsf1-regulated genes involved in protein folding, proteolysis, and reaction to heat (Albert et al., 2019). In our sc profiles, Sd cells exhibit increased expression of RASTR upregulated genes and constitutive targets of Hsf1 (Figure 4A and B), and we have previously shown (Figure 3E–F) that Sd cells exhibit decreased expression of ribosome processing and ribosome protein-encoding genes. This suggests that both the Rd and Sd states may contribute to a RASTR-like response in the days following FCZ exposure, and based on the RASTR response, that cells in the Rd state may transition to an Sd state.

Figure 4
Download asset Open asset
The Uniform Manifold Approximation and Projection (UMAP) embedding from Figure 3A, which contains fluconazole (FCZ) day 2 and day 3 cells.

Here, each Ribo-dominant (Rd) and Stress-dominant (Sd) cell is colored by its signature score for: (A) C. albicans orthologs of genes upregulated in ribosome assembly stress response (RASTR) and (B) constitutive targets of Hsf1.

Table 1
(A) Orthologous genes in C. albicans associated with ribosome assembly stress response (RASTR) in S. cerevisiae.

(B) Constitutive Hsf1 target genes.

RASTR signature (Albert et al., 2019)
DOWNRibosome processingASC1, RPL10, RPL10A, RPL11, RPL12, RPL13, RPL14, RPL16A, RPL17B, RPL18, RPL19A, RPL2, RPL20B, RPL21A, RPL23A, RPL24A, RPL25, RPL27A, RPL28, RPL30, RPL32, RPL35, RPL37B, RPL38, RPL39, RPL42, RPL5, RPL6, RPL8B, RPL9B, RPP0, RPP1B, RPS1, RPS10, RPS12, RPS13, RPS14B, RPS15, RPS16A, RPS17B, RPS18, RPS19A, RPS20, RPS21, RPS21B, RPS22A, RPS23A, RPS25B, RPS26A, RPS27, RPS27A, RPS3, RPS5, RPS6A, RPS7A, RPS8A, RPS9B, UBI3, YST1
UPProtein folding, response to heat, proteolysisUBI4, RPN4, PIN3, STF2, KAR2, MSI3, HSP60, HSP90, HSP70, SIS1, SSA2, HSP104, HSP78, STI1
Glucose and pyruvate metabolic processCYP1, PMA1, GLK1, TDH3, CDC19, PGK1
Unknown functionMBF1, KRE30, YDJ1, YBN5, ACT1, UBC4
Hsf1 constitutive targets (Leach et al., 2016)
ACE2, ADAEC, AHA1, ALO1, ALS1, ALS3, ALS4, ASR1, BOI2, BUL1, CCP1, CDC37, CDC48, CPR6, CRD2, CTF1, CYC1, CYP1, ERG2, GAP1, GIT3, GLX3, GOR1, GPX2, GPX3, GRP2, HCH1, HSF1, HSP104, HSP21, HSP60, HSP70, HSP78, HSP90, IFA14, ITR1, KAR2, MDJ1, MGE1, MIA40, MNN24, MSI3, PDC11, PGA56, PGA62, POR1, ROB1, RPM2, RPN4, RPS27A, SBA1, SBP1, SGT2, SIS1, SOK1, SSA2, SSC1, STI1, TRX1, TSA1, TSA1B, YDJ1, YWP1, ZCF35, ZPR1, ZWF1

Ifh1 is also a key component of RASTR control and signaling (Albert et al., 2019), and the vast majority of the targets for this transcription factor are RP genes (N=144 genes of which 41 are putative or uncharacterized genes; 64 RP) (Wade et al., 2004). Of these, 103 are detected at significant levels in our sc profiles. As shown previously, the RP genes are highest in Rd and lowly expressed in Sd.

The Rd and Sd subpopulations persist at 6 days post-FCZ treatment

Next, we explored whether the Rd and Sd subpopulations would persist once they were no longer exposed to the drug. On day 3 the FCZ-exposed cultures were transferred to YPD without drug and sc profiling was conducted 3 days later (FCZ day 6, Figure 1A, Materials and methods, ‘Strains, media, and drug treatment’). Figure 5A shows the UMAP embedding of the sc transcriptional profiles, using only the FCZ day 3 and FCZ day 6 survivor populations.

Figure 5
Download asset Open asset
Both the Rd and Sd states persist at six days post-fluconazol treatment.

(A) Uniform Manifold Approximation and Projection (UMAP) embedding of fluconazole (FCZ)-treated cells at days 3 and 6 (after resuspension in fresh media at day 3). Leiden clustering identified four clusters. Clusters are related along the trajectory (black) from slingshot analysis. The topology suggests that FCZ Ribo-dominant (Rd) cells at day 6 evolve from FCZ Rd cells at day 3. This relationship is also true for FCZ Stress-dominant (Sd) cells. (B) The heatmap depicts cluster stability analysis for the four clusters. Purple ticks depict the number of times each cell was assigned to each of the four clusters. Cluster stability analysis failed to strongly separate clusters 1 and 2. However, both of these clusters correspond to Rd cells. (C–F, H) The UMAP embedding of cells from (A) was annotated with color bars in C–F, H to depict: (C, F) density of cells from day 3 in FCZ and (F) from day 6 (after 3 days in YPD alone); or (D, E, H) signature scores for the (D) ribosomal protein (RP), (E) heat-shock stress, and (H) alcohol dehydrogenase gene signatures. (G) A heatmap depicting the average level of activation (VISION z-scores, color bar) of different signatures across clusters presented in (A), analogous in methodology to Figure 2C. (I) Growth curves for untreated (UT) and FCZ day 3 survivor resuspended in fresh media as measured by OD600 across days.

Leiden clustering was used to identify clusters of similarly behaving cells. Cluster stability analysis however did not find strong support separating clusters 1 and 2 (panel B), meaning that they could be collapsed together. FCZ day 3 cells are primarily located on the right in clusters 2 and 4; FCZ day 6 cells are primarily located on the left in clusters 1 and 3 (panels C and F). We detected the Rd response, predominantly in clusters 1 and 2. Specifically, these two clusters have RP and low HSP scores (Figure 5D and E). The Sd response is predominantly in clusters 3 and 4. Note that both the Rd and Sd subpopulation are detected in FCZ day 6 cells (cluster 1 and cluster 3, respectively). This suggests that the bifurcation is still present 3 days post-drug treatment after receiving fresh growth media.

Furthermore, because clusters 1 and 2 represent cells with only minor differences in the underlying expression patterns (because the two clusters are not stable), and because day 3 FCZ cells are primarily in cluster 2 and day 6 FCZ cells are primarily in cluster 1, we can conclude that day 3 and day 6 FCZ Rd cells have strong similarities. Given the strong resemblance in expression patterns between the Rd cells at day 3 and day 6, it is likely that the day 6 cells are ‘descendants’ of day 3 cells. Moreover, these Rd cells largely remain in the same state regardless of the increased time since FCZ exposure, suggesting Rd cells remain viable or largely unchanging at the transcriptional level.

Clusters 3 and 4 harbor day 6 FCZ and day 3 FCZ Sd cells, respectively. There are significantly differentially expressed genes (DEGs) between these two groups for example, the alcohol dehydrogenase signature is most highly expressed at day 6 (Figure 5G and H). Heat-shock stress is slightly reduced in day 6 compared to day 3 (Figure 5G). We also observe a substantial change in the relative proportions of the Rd and Sd populations: at day 3 we observe 70% Rd and 30% Sd but at day 6 we observe an almost even ratio (56% Rd and 44% Sd) (Figure 5C versus Figure 5F). Upon resuspension in fresh media, OD600 measurements confirmed that cells exposed to FCZ for 3 days and then resuspended in fresh media and allowed to grow for 3 days without drug, proliferate more slowly, and reach lower biomass levels than UT cells propagated without drug for 3 days (Figure 5I). Together this suggests that Sd cells, unlike Rd cells, are significantly changing over time since the FCZ exposure (although they do retain many of the molecular properties of day 3 cells) and may have either a higher proliferative capacity or survival rate relative to Rd cells.

We also used computational trajectory analysis to investigate the relationship between these clusters (Figure 5A, black tree; Materials and methods, ‘Cell clustering, trajectory, and signature analyses’, slingshot analysis). Slingshot analysis partially orders clusters of cells based on their transcriptional profiles; the resultant order is suggestive of how cells evolve their expression patterns over time in a manner analogous to an animation. Overall, these findings suggest that there is more evidence that day 6 FCZ Sd cells are ‘descendants’ of day 3 FCZ Sd cells compared to, for example, day 6 Sd cells ‘descending’ from day 3 Rd cells.

Discussion

Sc approaches can facilitate the characterization of phenotypic heterogeneity in microbial populations by enabling the detailed examination of individual microbial cells, revealing the full range of diversity in their gene expression. Here, we characterized two distinct subpopulations of cells with different expression profiles – the Rd and Sd states – among isogenic cells treated with FCZ, the most widely used antifungal drug. This result supports the observation that antifungal tolerance is due to a subpopulation of cells that grows in the presence of FCZ, while other cells in the population do not.

Despite the technical challenges associated with sc profiling of yeast cells, we identified heterogeneity in the gene expression of UT prototrophic SC5314 C. albicans cells that reflects the cell cycle phase and associated metabolic and stress responses. Specifically, glycolysis and several heat-shock genes were more highly expressed during M phase, and oxidative stress genes were more highly expressed in a distinct subset of cells that appear to be from G1/S phase. The existence of these distinct subpopulations is consistent with observations in the literature (Chiu et al., 2011; Finkel and Hwang, 2009; Hossain et al., 2021; Senn et al., 2012). The magnitude of stress signatures is much lower than in cells exposed to antifungal drugs. We also validated the distinct subpopulations by fluorescence microscopy with tagged genes differentially expressed in different cell cycle stages. Thus, the sc transcriptional profiling detects cell-to-cell heterogeneity within UT C. albicans cultures.

We then examined the transcriptional profiles of individual drug-treated cells. During FCZ exposure, distinct subpopulations of C. albicans appeared: on day 2, FCZ-treated cells were either in the Rd or Sd state with a higher frequency of Rd cells. By day 3, the Sd subpopulations become more balanced suggesting either the Sd subpopulation is more proliferative, or perhaps that cooperation between the two subpopulations is required, for example, for the exchange of metabolites. Before and after treatment, the vast majority of cells are in the white morphology so C. albicans morphological state does not explain this Rd/Sd population split.

The Rd response is characterized by high expression of ribosomal protein (RP) and ribosomal RNA (rRNA) processing genes. An imbalance of RP and rRNA can lead to proteotoxic stress and RP aggregation in S. cerevisiae (Tye et al., 2019), and might do so as well in C. albicans. The Sd state involves the heat-shock transcription factor Hsf1, which upregulates HSP90, HSP70, and other genes in the UPR, which have been shown to suppress RP expression in S. cerevisiae to protect against toxic protein build-up, and restores cell growth (Albert et al., 2019). The molecular profile of Sd cells is consistent with this transition exhibiting high expression of RASTR-induced genes, for example, Hsf1 constitutive target genes including HSP70, HSP90, and disaggregases HSP21, -60, -78, and -104, along with low expression of RP genes and rRNAs. Importantly, Sd cells also highly express several genes associated with drug tolerance (Delarze et al., 2020; Garnaud et al., 2018; Liu et al., 2015; Mayer et al., 2013; Nagao et al., 2012; Rosenberg et al., 2018) as well as genes involved in filamentous growth, adhesion, and biofilm formation. This is consistent with a potential role of Hsf1 (Hahn et al., 2004; Leach et al., 2016; Leach et al., 2012; Veri et al., 2018) as a trigger of RASTR.

Overall, the molecular profiles of the Rd and Sd states suggest a model in which some fraction of the cells exposed to FCZ initially enter the Rd response and trigger the RASTR response (Albert et al., 2019) which in turn causes the cells to transit into the Sd state. This could suggest that Sd cells may represent a C. albicans subpopulation that have mounted a successful drug tolerance stress response. Alternatively, it could be the case that both subpopulations are necessary for cell survival.

We also treated cells with RAPA (analysis in Appendix 2). RAPA in isolation is considered to have weak antifungal properties, but drug tolerance (but not drug resistance) is ablated when used in combination with FCZ (Rosenberg et al., 2018). We did not identify any RAPA-treated cells in the Rd state. This is perhaps not surprising given that RAPA targets Tor1 and shuts down ribosomal biogenesis, an important component of the Rd state. In fact, RAPA-exposed cells also had the weakest RP expression. Of note, RAPA specifically triggers the expression of genes involved in hyperosmotic stress. Previous work established that the Tor1 pathway controls phase separation and the biophysical properties of the cytoplasm by tuning macromolecular crowding (Delarue et al., 2018). Hyperosmotic stress response to RAPA exposure might reflect disruption of the biophysical properties of the cytoplasm, altering the osmotic balance in cells.

Whereas nearly every RAPA-treated cell had an Sd transcriptional profile, nearly every CSP-treated cell after 2 days exposure had an Rd-like transcriptional profile (analysis in Appendix 2). The echinocandin CSP is fungicidal with a mechanism of action distinct from FCZ – it inhibits β-1,3 glucan synthesis and disrupts the fungal cell wall. We did not find any CSP-treated cells in the Sd subpopulation by day 3. Although we used low concentrations of CSP to avoid cell aggregation and therefore it is possible that alternative CSP concentrations could induce Sd cells, it appears that survivors of the fungicidal drug are transcriptionally far more homogeneous at early time points (days 2–3).

After resuspension in fresh media and 3 days of growth (day 6), both the Sd and Rd subpopulations persist. In fact, now Sd outnumbers the Rd subpopulation. Their persistence hints that these surviving cells have some ‘memory’ of treatment with FCZ. In the case of Rd cells, there are few transcriptional changes between days 3 and 6. For Sd, there is a downregulation of the heat-shock response and an increase in the alcohol dehydrogenase pathway, suggestive of cellular proliferation in the fresh medium. Sd cells at day 6 still however largely retain differential expression of genes, pathways, and proteins identified in Sd cells at day 3.

Conclusions

The use of a nanoliter droplet-based assay adapted for fungal cells enabled a detailed sc analysis of C. albicans in the absence and presence of antifungal drugs over a period of days. The assay was cost-effective and encountered minimal issues or failed runs. Studies of thousands of individual cells enriched our understanding of community structure and population heterogeneity. Specifically, this study refines published bulk transcriptome studies by differentiating between genes, pathways, and responses that are expressed universally in all cells, and those that are restricted to specific subpopulations. Here, by examining cellular trajectories across populations in a high-throughput manner, we obtained new insights into time-sensitive processes, such as the emergence of drug tolerance during 2–3 days of drug exposure and the existence of cells with Rd and Sd states, which reflect two different cellular states within the RASTR response described as a molecular signature of an S. cerevisiae stress response. If the Rd response proves generalizable to other drug treatments, disrupting the Rd to Sd tolerance transition or the Sd response itself could represent an innovative therapeutic approach for other antifungal treatments.

Materials and methods

Strains, media, and drug treatment

C. albicans cultures for sc-RNA-, bulk RNA- and DNA-sequencing

Request a detailed protocol

C. albicans SC5314 cells were streaked out from glycerol stocks in –80°C on YPD agar plates (2% D-glucose, 2% peptone, 1% yeast extract, 0.01% uridine, 2% agar) and incubated at 30°C for 48 hr. Afterward, a single colony of cells was transferred into YPD liquid media (2% D-glucose, 2% peptone, 1% yeast extract, 0.01% uridine) and incubated at 30°C for 12–16 hr.

Preparation for sc profiling via DROP-seq

Request a detailed protocol

For UT cells, a single colony of SC5314 cells were transferred to 5 ml YPD and grown overnight to yield ~108 cells/ml. Cultures were then diluted to an OD600 of 0.1 in fresh 50 ml YPD liquid and incubated at 30°C in a shaker incubator. Cells were collected when OD reached 0.5–0.9 in order to maximize the number of cells in mid-log phase. Cells were pelleted by centrifugation, 1 ml of RNAlater was added (Sigma # R0901), the suspended cells were incubated for 10 min at room temperature, and the resulting culture was frozen at –20°C for later use.

For treated cells, we performed the following approach to ensure a sufficient number of survivor cells for profiling across different drugs and time points: cells were pelleted and resuspended in 1 ml of YPD, 250 μl of this suspension was combined with 15 ml of fresh YPD and incubated at 30°C with shaking until the culture reached an OD600 of 0.5–0.6. Finally, ~108 of these cells were seeded into 10 ml of YPD. Each suspension was then subjected to drug treatment.

FCZ (Sigma #F8929) was used at 1 µg/ml, which is from 1× to 2× the dosage relative to the MIC50 for SC5314 in YPD (Figure 3—figure supplement 1). CSP (Sigma #SML0425), a compound that interrupts cell wall biosynthesis (Cappelletty and Eiselstein-McKitrick, 2007; Stevens et al., 2004; Yang et al., 2017), was used at 1 ng/ml, which is well below the reported MIC50 and was chosen to ensure a sufficient number of non-aggregated survivors to generate sc profiles. For RAPA, a subinhibitory concentration of 0.5 ng/ml was chosen based on previous studies that established such levels are sufficient to generate a fungicidal synergistic effect when given concomitantly with FCZ (Rosenberg et al., 2018). Each drug was added to the individual cultures and incubated at 30°C for 48 or 72 hr. For the day 6 population, day 3 survivors were washed twice, and resuspended in 10 ml of fresh YPD.

At each time point, cultures were collected and strained (pluriStrainer 20 µm, pluriSelect) before placement in fresh tubes. Straining was done in order to minimize the likelihood of clogging in the microfluidic due to rare but large hyphae and pseudohyphae morphologies. We observed that germ tubes up to four times the length of the mother cell can still be processed for DROP-seq analysis. Such cells are well within the hyphal transcriptional profile (Nantel et al., 2002), suggesting that our results may contain some profiles of hyphae and pseudohyphae cells. After filtering, the vast majority of cells were in the yeast white morphology with less than 0.2% of cells in a filamentous morphology (hyphae or pseudohyphae) after manual counting ~100 microscopy images with an average of ~50 cells per slide for each such population. All cultures yielded a sufficient population of survivors for downstream sc transcriptional profiling, bulk transcriptional profiling, bulk DNA genomic profiling, and/or microscopy. Cultures were washed with 1 ml of RNAlater twice. Cells were then resuspended in 1 ml RNAlater and incubated at room temperature for 10 min before storage at –20°C until sc profiling with DROP-seq.

Cultures for OD600 analyses

Request a detailed protocol

Three cultures were inoculated with single cultures and incubated overnight in YPD at 30°C. The following morning, cultures were diluted to an OD600 of ~0.5 in either fresh YPD or YPD with 1 µg/ml FCZ (same concentration used for the DROP-seq experiment). OD600 was measured immediately without dilution, then every 24 hr via a 1:100 dilution for a total of 72 hr (Figure 3—figure supplement 1A). To investigate the growth of day 3 FCZ survivors, three cultures previously grown for 72 hr in either YPD or YPD with 1 ug/ml FCZ were diluted to an OD600 of ~0.5 in fresh YPD without drug. OD600 was measured immediately without dilution and then every 24 hr via a 1:100 dilution for a total of 72 hr (Figure 5I).

Spheroplasts

Request a detailed protocol

The C. albicans setting required an optimized protocol for the removal of the cell wall and to induce stable spheroplasts for sc profiling. Toward this end, we experimented with different concentrations of zymolyase (0.1, 0.2, and 0.4 U zymolyase (BioShop # ZYM002) with 107 cells in 100 µl of sorbitol 1M) at different time points (incubated at 37°C for 10, 20, 30 min) before processing with the DROP-seq. To compare cultures grown under different conditions, cells were stained with calcofluor white and imaged using a Leica DM6000 microscope. We concluded that concentrations in the range 0.15–0.25 U after 25 min are able to induce spheroplasts that remain sufficiently stable for processing with our DROP-seq.

Sample preparation for sc profiling

Request a detailed protocol

At the time of DROP-seq profiling, an aliquot of 107 (OD600=0.68) cells was separated from the culture in Materials and methods, ‘Strains, media, and drug treatment’, and washed three times with sorbitol 1M solution. The cells were then resuspended in 100 μl sorbitol 1M+0.25 U zymolyase and incubated at 37°C for 25 min (as per our findings in Materials and methods, ‘Spheroplasts’). Next, the cells were pelleted and resuspended again in 0.5 ml of cold and fresh RNAlater for 5 min. Now, the cells were washed (centrifuged and pelleted) with 1 ml of washing buffer (1 M sorbitol, 10 mM TRIS pH 8, 100 µg/ml BSA) three times. Finally, 106 cells (OD600=0.08) were resuspended in 1.2 ml of the washing buffer. This cell suspension was then used as input to the DROP-seq device.

Sample preparation generally follows the protocol given by Macosko et al., 2015, with some exceptions. Whereas Macosko et al. recommend a ratio of 100K mammalian cells to 120K beads for DROP-seq, we found that a ratio of 1M cells for 120K beads generated a sufficient yield of cDNA as per the Tapestation (Agilent Inc) device. Jackson et al., 2020, use 5M cells as input to the Chromium (10X Inc) system. Furthermore, whereas Macosko et al. use 1 ml of lysis buffer, we use 1.2 ml. Instead of 13 PCR cycles, we use 17 (Jackson et al. uses 10 cycles). Samples were sequenced using a NEXT-seq 500 (Illumina Inc) following the standard Macosko et al. protocol set to yield an estimated 200 million reads/sample.

Quality control, basic processing, and normalization of the sc profiles

Request a detailed protocol

All computations were performed using Python version 3.9.6 (van Rossum and Drake, 2009) or R version 4.0.4 (R Development Core Team, 2021). Gene abundances were estimated from raw sequencing data using the end-to-end pipeline alevin-fry (He et al., 2021) which performs UMI deduplication and reduces the number of discarded (multimapped) reads. The pipeline utilizes a reference index covering the spliced transcriptome extracted from the latest version of C. albicans strain SC5314_A22 (haplotype A, version 22; GCF_000182965.3). The Unspliced-Spliced-Ambiguous (USA) mode was used to separately keep track of the types of transcripts from which UMIs are sampled. A gene-by-cell matrix for each sample was obtained by summing reads labeled as either ‘spliced’ or ‘ambiguous’ by alevin-fry.

We started by filtering genes from downstream analyses with a zero sum count (the count across all cells), as were cells with less than 5, or more than 2000 transcript counts (Figure 1—figure supplement 1A). Droplets almost always capture ambient ‘free floating’ RNA that is present in the suspension. Therefore, even when a microparticle is not captured alongside a cell in a droplet, sequencing still yields reads originating from this ambient RNA. However, the number of such reads is far lower than a droplet with a successful cell capture. Cells were included in their analysis if their profiles significantly deviated from levels indicative of ambient RNA in the suspension using EmptyDrops (FDR <0.01) (Lun et al., 2018). SCANPY (Wolf et al., 2018), a python-based toolkit, and SingleCellExperiment (Amezquita et al., 2020), a R/Bioconductor package, were used for data quality control, filtering of genes and downstream visualization. We removed genes (n=869) that were expressed in less than 20 cells (Figure 1—figure supplement 1C).

We then use scVI (Gayoso et al., 2021; Lopez et al., 2018) version 0.12.2, a Bayesian deep neural network architecture which implements a probabilistic model of mRNA capture and uses a variational autoencoder to estimate priors across batches and conditions. Models were built using default parameters for different grouping of samples: (i) UT cells, (ii) UT cells and CSP, RAPA, FCZ treated at day 2 (and 3), and (iii) FCZ-treated cells at days 3 and 6. All models were adjusted for batch and library size. We trained scVI’s variational autoencoder and stored the latent representation for visualization and downstream analyses. We reduced the inferred latent spaces to two dimensions via the UMAP tool using the implementation of umap-learn (McInnes et al., 2020) in SCANPY (Wolf et al., 2018) (min_dist = 0.3). When analyzing the level of expression of individual genes, missing or dropout values were first imputed using MAGIC (van Dijk et al., 2018).

Bulk transcriptomics

Request a detailed protocol

Total RNA was extracted from FCZ-treated cells at day 2 post-exposure, which were grown according to Materials and methods, ‘Preparation for sc profiling via DROP-seq:’, using the QIAGEN RNeasy mini kit protocol. RNA quality and quantity were determined using a Bioanalyzer (Agilent Inc). Paired-end read sequencing (2×50 bp) was carried out on a NextSeq500 sequencer (0.5 Flowcell High Output; Illumina Inc). Raw reads were pre-processed with the sequence-grooming tool cutadapt (Martin, 2011) version 0.4.1 with quality trimming and filtering parameters: --phred33 --length 36 -–2colour 20 --stringency 1 -e 0.1. Each read pair was mapped against C. albicans strain SC5314_A22 (haplotype A, version 22; GCF_000182965.3) downloaded from the NCBI using STAR (Dobin et al., 2013) version 2.7.9a with the following filtering parameters: --outSAMmultNmax 1 --outSAMunmapped Within --outSAMstrandField intronMotif. We obtained ~13 million reads of which 88% were uniquely mapped along the genome. The read alignments and C. albicans genome annotation strain SC5314_A22 (haplotype A, version 22; GCF_000182965.3) were provided as input to featureCounts() from the Rsubread package (Liao et al., 2019) version 2.4.3 to estimate gene abundances. The following parameters were used: isPairedEnd = TRUE, countReadPairs = TRUE, requireBothEndsMapped = TRUE, checkFragLength = FALSE, countChimericFragments = FALSE, countMultiMappingReads = TRUE, fraction = TRUE.

Construction of pseudo-bulk profiles

Request a detailed protocol

Throughout the manuscript, pseudo-bulk profiles refer to transcriptional profiles that are derived from the sc reads by ignoring barcodes. This pipeline is depicted in Figure 1—figure supplement 2A. By ignoring the R1 (left) read of the sc profile that contains the cellular barcode, we are effectively performing ‘bulk’ RNA-sequencing using only the R2 (right) read that aligns to a transcript in the sample. This collapses all cells to a single profile. We compared two different techniques to compute pseudo-bulk profiles, or we can first use the barcoded reads to partition cells into classes (e.g., Rd versus Sd) and then form a pseudo-bulk profile specific to each of the classes. The unfiltered pseudo-bulk data is derived from counting raw reads aligned to the reference genome using the STAR tool (Dobin et al., 2013). The filtered pseudo-bulk dataset is obtained by first applying our sc pipeline (alevin-fry followed by EmptyDrops) and then summing across all cells. The first is closer in spirit to true bulk (single read) profiling, while the second approach filters reads, cells, and genes in the same manner as sc analyses and therefore represents a middle point between bulk and sc profiling. A comparison of filtered pseudo-bulk FCZ profiles at day 2 and 3 datasets indicated that the assay is robustly quantifying the expression of genes across different batches (Figure 1—figure supplement 2A).

A comparison of bulk profiles versus unfiltered pseudo-bulk for the FCZ profiles at days 2 and 3 indicated that both methods identified all but 297 of the same genes. The missed genes tended to be expressed at low levels in the bulk profiles (Figure 1—figure supplement 2B). Moreover, day 2 and day 3 pseudo-bulk profiles were significantly correlated with ‘true’ bulk RNA-sequencing profiles (R ranges from 0.67 to 0.74; Figure 1—figure supplement 2C). Here ‘true’ bulk RNA-sequencing profiles were generated as described in Materials and methods, ‘Bulk transcriptomics’.

Whole-genome DNA-sequencing

Request a detailed protocol

C. albicans populations were grown as described in Materials and methods, ‘Strains, media, and drug treatment’, although we did not apply a cell filtration step to remove filamentous cells. Preparation of genomic DNA used the MasterPure Yeast DNA Purification Kit (Lucigen # MPY80200) with a NextSeq500 – 1 flowcell mid output (130M fragments), 150 cycles pair-end reads (maximum 2×80 nt), yielding on average 23.9 million reads per sample (1 UT, FCZ at days 2, 3, 6, and 12). This gives an expected sequencing depth of 224 since the size of the C. albicans genome is ~16 Mb.

Raw reads were pre-processed with the sequence-grooming tool cutadapt (Martin, 2011) version 0.4.1 with quality trimming and filtering parameters: --phred33 --length 36 -q 5 --stringency 1 -e 0.1. Analyses followed a previous study by Ford et al., 2015, which also examined C. albicans complete genomes, although we used a more recent version A22 of the C. albicans SC5314 haplotype A genome for read mapping (downloaded from the Candida Genome Database; https://www.candidagenome.org/). Briefly, reads were mapped using the BWA alignment tool (Li and Durbin, 2009) and RealignerTargetCreator and IndelRealigner from GATK (McKenna et al., 2010) were used for re-alignment. SNPs were detecting using Unified Genotyper (GATK version 1.4.14) using the same filtration criteria as reported in Ford et al. Determination of copy number and loss-of-heterozygosity was performed using GATK with the strategy reported in Ford et al.

Restricting attention to the UT samples, our observed sequencing depth was just under 200 and we identified approximately the same number of single nucleotide polymorphisms (N=3304) and the same number of insertions/deletions (181/255 resp.) as a previous whole-genome sequencing effort (Cuomo et al., 2019) using their bioinformatic pipeline. Although mutations were observed in some reads at some genomic loci, the samples were not enrichment for mutations that occurred more often than the rate of sequencing error which is 10–3 after correcting for multiple testing. This error rate is in line with estimates of the spontaneous mutation rate for C. albicans (Ene et al., 2018), together suggesting that the population is near isogenic without any significantly large subclones. Given that the error rate for copy number variants (e.g. loss, amplification) is lower than polymorphisms (Ene et al., 2018) and the duration of cell expansion before drug exposure (<2 days), it would be unlikely that spontaneous mutations explain the degree of heterogeneity that was observed 2–3 days post-drug exposure.

Cell clustering, trajectory, and signature analyses

Request a detailed protocol

To identify subpopulations of cells with similar gene expression patterns in an unbiased, unsupervised manner, we applied Leiden clustering (Traag et al., 2019) on the latent space generated by scVI (resolution of 0.5 for the analysis combining UT, FCZ day 2 and 3, CSP, and RAPA-treated cells and resolution of 0.4 for FCZ day 3–6 analyses). To test the robustness of the clustering process, we repeated the clustering process 100 times, each time using a random subset of 95% of the cells. With each such random set, we repeated the scVI model building, followed by Leiden clustering to identify clusters. The mapping from the original clusters to the new clusters was then established based on the maximum overlap in cell membership between the original and newly formed clusters. Finally, we counted the frequency that each cell was assigned to each cluster over all the iterations.

The lineage trajectory of day 3 and 6 FCZ-treated cells was performed using the R/Bioconductor slingshot package version 1.8.0 with default parameters (Street et al., 2018).

To investigate the key sources of variability across cell subpopulations, we curated 43 gene signatures related to microbial phenotypic diversity and drug tolerance including the cell cycle, TCA cycle, specific and general stress responses, metabolic pathways, amino acid biosynthesis, efflux pumps, and specific drug responses among others (Azadmanesh et al., 2017; Berman, 2006; Côte et al., 2009; Enjalbert et al., 2006; Enjalbert et al., 2003; Gasch et al., 2000; Hardwick et al., 1999; Hinnebusch, 2005; Jackson et al., 2020; Natarajan et al., 2001; O’Duibhir et al., 2014; Pais et al., 2019; Sanglard, 2016; Sanglard et al., 2003; Tsai et al., 2019; Yang et al., 2017; Figure 2—source data 1). In some cases, gene signatures from the literature were derived in other organisms and required orthology mappings to C. albicans.

Briefly, signatures of cell cycle phases were identified as transcriptional expression patterns in synchronous C. albicans populations (Berman, 2006; Côte et al., 2009) or were expert-curated and well-established cell cycle genes found in distinct clusters of S. cerevisiae sc profiles (Jackson et al., 2020). Signatures specific to certain stresses were found by transcriptional profiling of C. albicans challenged by temperature, osmotic and oxidative stress under conditions that permitted >60% cell survival (Enjalbert et al., 2003) or in C. glabrata (Pais et al., 2019). We also curated more general non-specific stress signatures identified in C. albicans (Enjalbert et al., 2006) or S. cerevisiae (Gasch, 2007; Gasch et al., 2017; Tsai et al., 2019). Previous studies established the existence of a ubiquitous environmental stress response (ESR) in S. cerevisiae which is deregulated in response to many different environmental perturbations (Gasch, 2007; Gasch et al., 2000). The ESR is divided into the induced (iESR) and repressed (rESR) subcomponents. The iESR is characterized by overexpression of heat-shock and oxidative stress genes in addition to genes involved in central carbohydrate metabolism and energy generation (Gasch, 2007; Gasch et al., 2000). Previous studies have noted a complex, intricate relationship between the ESR and cell cycle phase (O’Duibhir et al., 2014; Regenberg et al., 2006). Some investigations suggested a smaller core ESR in C. albicans (Enjalbert et al., 2006) which may have evolved due to the unique host environment with the need to grow with different substrates (Brown et al., 2014). Curated metabolic signatures include lowly expressed glycolytic genes and highly expressed TCA cycle genes identified during diauxic shift or RAPA treatment in yeast (DeRisi et al., 1997; Hardwick et al., 1999) as well as GCN4-driven amino acid biosynthesis, another well-defined metabolic signature conserved between S. cerevisiae and C. albicans (Hinnebusch, 2005; Natarajan et al., 2001).

The signature analyses start with the VISION tool which estimates a signature score for every cell (DeTomaso et al., 2019) using the batch-adjusted normalized counts returned by scVI’smodel. The distribution of individual scores for cells classified in each cluster can be depicted using the empirical cumulative distribution function. These distributions were further compared using the Kolmogorov-Smirnov test. We then used the median to summarize signatures scores of cells within each cluster and selected signatures which were the most variable across clusters (sd >0.05 for the analysis combining UT, FCZ- day 2 and 3, CSP-, and RAPA-treated cells or combining FCZ-treated cells at days 3 and 6; ). Heatmaps were used to depict the median scores (z-score, color bar) of the selected signatures for each cluster.

DGE analysis

Sc differential gene expression (DGE)

Request a detailed protocol

scVI’s model allows us to approximate the posterior probability of the alternative hypotheses (genes are different) and that of the null hypotheses through repeated sampling from the variational distribution, thus obtaining a low variance estimate of their ratio (i.e., Bayes factor). We used this approach to identify genes differentially expressed in the combined and individual comet clusters compared to the other clusters (Appendix 2—figure 2A, B).

Pseudo-bulk DGE and gene ontology (GO) enrichment analysis

Request a detailed protocol

We also conducted pseudo-bulk differential gene expression analyses which allow for a dramatic reduction in the number of zeros in the data by aggregating cells within each replicate. This approach has been found to achieve the highest fidelity to the experimental ground truths significantly reducing the risk of false discoveries (Squair et al., 2021). As described in Materials and methods, ‘Construction of pseudo-bulk profiles’, filtered pseudo-bulk profiles were obtained by summing counts of selected cells within each replicate. In order to use DESeq2 (Love et al., 2014), a standard R/Biocoductor package used for differential analysis of count data, we require at least two replicates within each group of comparison. If only one replicate was available (e.g., FCZ at day 6), we partitioned the selected cells of a single replicate into two groups randomly before forming pseudo-bulk profiles. We then used DESeq2 default parameters to perform DE analysis and selected genes with Benjamini-Hochberg FDR <0.1 (Benjamini and Hochberg, 1995). Models were adjusted for batch affects where relevant.

Finally, we identified enrichment of biological processes in lists of significantly over- and under-expressed genes using the R/Bioconductor ViSEAGO package (Brionne et al., 2019). ViSEAGO includes all algorithms developed in the R/Bioconductor topGO package including the weight01 fisher test that takes into account the topology of the GO graph (Alexa et al., 2006). Biological processes with weight01 p-value <0.01 were defined as significantly enriched in the gene list.

Cell imaging

Request a detailed protocol

To validate the subpopulation structure identified by the sc transcriptomics, we choose markers representative of distinct clusters. Cells were transformed with GFP and RFP fusion constructs for HSP70 and TTR1 marker genes respectively using a CRISPR/Cas9 protocol (Min et al., 2016) with primers described in Figure 2—source data 2. Strain SN76(his1Δ/his1Δ, arg4Δ/arg4Δ, ura3Δ/ura3Δ) was chosen for gene tagging since it is a derivative strain of SC5314 but with multiple auxotrophic markers. These markers allow for convenient selection of successfully transformed cells.

Benchling (https://benchling.com) was used to design the sgRNAs. We followed the CRISPR/Cas9 protocol with the plasmid pV1093 from Min et al., 2016. This includes two PCRs to fuse the SNR52 promoter to the sgRNA scaffold and terminator. The third PCR amplifies the final sgRNA cassettes. Two different plasmids pENO1-iRFP-NATr (Plasmid #129731, Addgene Inc) and pFA-GFP-HIS1 were used to design the repair segments. The construction of the Cas9 cassette proceeded as per Min et al. Amplification of the Cas9 cassette with PCR used the following schedule: 98°C for 3 min, 98°C for 30 s, 63°C for 30 s, 72°C for 5 min and 30 s. Steps 2–4 have been repeated for 34 rounds followed by 72°C for 10 min and finally the reaction finished in 4°C. The repair DNA must be amplified with the designed primers (Figure 2—source data 2) in 8–12 PCR tubes with 0.1 μl plasmid (500 ng/ml), 2.5 μl forward primer, 2.5 μl reverse primer, 1 μl 10 mM dNTP, 33.65 μl nuclease free water, 10 μl 5× HF PCR buffer, and 0.25 μl phusion polymerase in each tube.

Cells that were successfully transformed were grown and harvested in a manner identical to that used for the sc experiments (Materials and methods, ‘Strains, media, and drug treatment’). At time of microscopy, cells were collected, washed with H2O, and transferred to minimum media to minimize the background noise from normal YPD media. Afterward, cells were mounted onto the uSlide and imaged with a Leica DM6000 microscope at 1000× (~50 images/time point; ~50 cells/image).

Appendix 1

Drug diffusion assays with FCZ-treated cells at days 1–6

Diffusion disk assays provide a means to measure tolerance and resistance of a fungal population to a drug (Gerstein et al., 2016). In our case, we use these assays to measure tolerance and resistance of C. albicans in response to FCZ. Toward this end, two SC5314 colonies were transferred from YPD agar plates to 5 ml YPD liquid media and incubated at 30°C overnight. Cells were washed the next day with PBS 1× twice and an aliquot of 106 cells/ml was transferred into an Eppendorf tube. 100 µl of this suspension was plated on YPD agar and spread evenly with a cell spreader. A FCZ disk (Bio-Rad #62802) was placed at the center of the plate and incubated at 30°C. This was repeated three times (replicates 1–3) for a total of two biological replicates (experiments 1–2).

Images were captured every 24 hr for 6 days (Appendix 1—figure 1). Appendix 1—figure 2 describes how the radius (RAD) and the fraction of growth (FoG) are estimated from these images to provide measures of drug resistance and tolerance, respectively. We use the diskimageR R package to estimate the intensity of 72 lines (in 5° increments around the circle) from the disk edge to the rim of the plate for each replicate across 6 days post-exposure (Appendix 1—table 1). A Kruskal-Wallis χ2 test was used to test whether there was a difference in the mean of either the FoG or RAD between day 1 and the remaining time points versus a null hypothesis of no difference.

Appendix 1—figure 1
Download asset Open asset
Fluconazole (FCZ) disk assay results across 6 days.

Two cultures (biological replicates E1 and E2) of C. albicans SC5314 were grown in YPD and spread on three plates (technical replicates R1, R2, and R3). The bright white spot at the center of the plate is the FCZ diffusion disk.

Appendix 1—figure 2
Download asset Open asset
Diffusion disk assay experiments with fluconazole.

(A) A C. albicans culture with fluconazole (FCZ) disk in the middle of plate. The red radial line represents 1 of 72 measurements taken every 5°. (B) Top is a blown-up and restricted image of the red line from (A) and bottom represents the pixel intensities from 0 (edge of disk) to 40 (edge of plate). (C) Gray curves are from the 72 measured radial lines. Black represents the average of these 72 measurements per mm from disk edge (x-axis). Light blue dot is the RAD20, defined as the point where there is a 20% reduction in growth. Middle blue dot is the RAD50 and dark blue represents the RAD80. The FoG20 is defined as the fraction of the area under the black line from 0 to the RAD20 (i.e., the amount of growth observed) divided by the total potential growth (delimited by dotted line). The FoG50 and FoG80 are defined analogously (adapted from Gerstein et al., 2016).

Appendix 1—table 1
Two biological (independent) C. albicans SC534 cultures and three technical replicates were subjected to fluconazole (FCZ) diffusion disk assay (DDA) experiments and averaged.

The table represents the output from diskimageR over 6 days of imaging. The increase in tolerance, which is measured by the fraction of growth (FoG), from day 1 to day 2 is significant as is the increase from day 2 to day 3 (both p<0.001, Kruskal-Wallis χ 2). There is a slight but statistically reduction in the RAD20, suggesting that there has been a small increase in resistance (p<0.05, Kruskal-Wallis χ2). However, we can see that the vast majority of the culture is not resistant.

DayFoGRAD20
10.1416.33
20.2114.17
30.3712.67
40.3713.00
50.3812.83
60.413.00

Appendix 2

Supplemental results

UT cells are found primarily in cluster 3-green and to a lesser extent in cluster 2-pink (Figure 3A–B; Appendix 2—figure 1A). CSP cells were mostly found in cluster 1-darkpink while RAPA cells were mostly found in cluster 2-lightpink (Figure 3A and B; Appendix 2—figure 1). The remaining 3% of cells (N=573) were outliers scattered across 19 distinct clusters (others cluster-darkblue; Figure 3A and B) that we termed ‘comets’. In this section, we will describe the biological characteristics for each of these clusters.

UT cells show expectedly higher glycolysis gene expression compared to antifungal-treated cells on days 2–3

UT cells have the highest expression of the glycolysis and carbohydrate reserve metabolism pathways (Figure 3C, green column; Appendix 2—figure 1). This is consistent with the fact that UT cells were collected during log-phase before glucose became a limiting factor, while all non-UT cells were harvested at least 2 days post-media renewal.

Exposure to FCZ: emergence of ‘comets’

In addition to the five clusters that comprise the vast majority of all profile cells, we observe 19 small clusters that are scattered across the UMAP embedding (Figure 3A). We term these 19 small clusters ‘comets’. Using the sc profiles, we investigated whether there is evidence that the comets are interesting biological phenomena versus stochastic noise from the sc profiling technique. Since the small size of the comets (comets contain 7–126 cells) creates significant technical/statistical challenges, our investigation and findings are built using several distinct analysis methods.

DEGs: all comets combined versus the five major clusters. All cells that belong to comets were combined into one ‘super comet’ and contrasted against the five major clusters to identify DEGs using two different approaches (Materials and methods, ‘DGE analysis’; Appendix 2—figure 2; Appendix 2—figure 2—source data 1). These analyses identified GCN4 (a transcriptional activator of the general amino acid response), and ILV2 and ILV5 (which encode for enzymes involved in branched-chain amino acid biosynthetic pathway) (Appendix 2—figure 2 ; Bayes factor >2.5 and proportion of non-zero value >0.2). The ILV pathway is known as a fungal amino acid starvation therapeutic target and ILV5 is recognized as a component in the C. albicans amino acid starvation response (Kingsbury and McCusker, 2010; Tripathi et al., 2002). It is regulated by the transcriptional activator GCN4 in both S. cerevisiae and C. albicans (Hinnebusch, 2005; Natarajan et al., 2001; Tripathi et al., 2002).

Pathway and signature analysis can sometimes detect more subtle differences in comparison to univariate DEG analysis, especially in cases such as ours where there are few cells available for the analysis. Figure 3C (column labeled ‘others’) summarizes our observations from this analysis, which identifies an enrichment of GCN4 target genes encoding amino acid biosynthetic enzymes, suggesting that comet cells are undergoing amino acid starvation. There are at least three reasons why C. albicans cells may be depleted of amino acids: (i) exposure of C. albicans to FCZ reduces the pool size of amino acid biosynthesis intermediates (Katragkou et al., 2016), (ii) the cells are exporting amino acids, which is consistent with the observed upregulation of efflux signatures in Figure 3C or (iii) the amino acid starvation response could be related to changes in the translational machinery when cells experience extreme stress. This is consistent with the observation of low expression of rRNAs in Figure 3C.

GO analysis of the DEGs also revealed an enrichment of genes involved in DNA replication regulation, actin cytoskeleton, plasma membrane organization, and mitotic ring assembly, all highly expressed in comet clusters (NHP6A, HHF22, HHF1, SBA1, NCE102, KIN2, RVS167, CDC12) (Appendix 2—figure 2—source data 1; 2; Appendix 2—figure 2). This is in line with the observation that FCZ affects bud formation and is associated with DNA replication or spindle pole body duplication (Harrison et al., 2014). In our prior research, we discovered that within 4–8 hr of drug exposure (at concentrations greater than the MIC), approximately 20% of mother-daughter cell pairs produced a single granddaughter cell, resulting in three cells (Harrison et al., 2014). Despite this, the mother and daughter nuclei continued to divide, generating four daughter nuclei. This process led to one cell containing twice the normal DNA content, and the resulting 4N cells frequently underwent abnormal divisions, giving rise to aneuploid offspring (Harrison et al., 2014). Using the sc data, we searched for aneuploidy evidence by calculating the percentage of reads assigned to genes in each chromosome. We found that cells in cluster 16 have a significantly higher proportion of reads assigned to genes in chromosome 2, which could indicate trisomy or local chromosomal amplification (Appendix 2—figure 2).

DEGs: each comet versus each main cluster. Although the number of cells per comet cluster is small (from 7 to 126 cells), individual molecular marker genes could be identified in some comets. Some of these markers have been reported previously to be involved in FCZ tolerance or resistance (Appendix 2—figure supplement 2B; Materials and methods, ‘DGE analysis’; top 5 with Bayes factor >3 and proportion of non-zero value >0.2). These included cell wall organization and adhesion (ADA2, TPK2, LMO1, KRE9, RDH54, DFI1) and chromosomal and cytoplasmic division (AXL2, HOS3, RDH54, RFA2 RIM15, NPL4, CDC5 BUB2).

In summary, although the comets each contain a small number of cells make it difficult to assess any statements with statistical rigor, comets generally overexpress GCN4 and related genes encoding amino acid biosynthetic enzymes, suggesting the cells are starved for amino acids. Pathway analysis is consistent with various known reasons for amino acid starvation. Specifically, efflux pumps are concomitantly overexpressed and translational machinery is underexpressed. Comets also consistently overexpress DNA replication regulation, actin cytoskeleton, plasma membrane organization, and mitotic ring assembly processes. We have previously established that FCZ exposure is associated with DNA replication or spindle pole body duplication, leading to aneuploidy. At least one comet has strong evidence of an aneuploidy along chromosome 2.

Subinhibitory dose of CSP triggers the Rd response

CSP, a first-line clinical treatment echinocandin, acts by inhibiting β-1,3 glucan synthesis and disrupting the fungal cell wall (McCormack and Perry, 2005). Although it is effective against most Candida species, a known mechanism of echinocandin resistance involves point mutations in ‘hot spot’ regions of FS-encoded subunits of a glucan synthase (Perlin, 2015; Pristov and Ghannoum, 2019). Tolerance to CSP increases with larger amounts of cell wall chitin, potentially due to aneuploidy and chromosome 5 rearrangement (Yang et al., 2017). Even though fungicidal CSP and fungistatic FCZ have different modes of action, they share similar tolerance mechanisms associated with Hsp90, the calcineurin pathway (Singh et al., 2009) and the pH-responsive RIM pathway (Garnaud et al., 2018). To treat a C. albicans population, we used a subinhibitory concentration of CSP, well below its MIC50 (1 ng/ml), as higher doses led to significant cell aggregation, rendering the cells unsuitable for sc profiling.

On day 2, the majority of cells were found in cluster 1-darkpink, where Rd FCZ cells also reside (86% of all CSP, Figure 3A and B; Appendix 2—figure 1). We examined the differences between CSP and FCZ cells in the Rd state and identified only one significantly differentially expressed gene (FDR <0.1), suggesting a highly similar Rd response for both antifungal treatments, at least at the selected dosage levels and day 2 time point. Unlike FCZ, no bifurcation was observed for CSP at day 2 (i.e., we did not observe any Sd cells), which could be partially due to our choice of a low dosage level, a factor previously shown to be crucial in determining the cellular response to CSP in C. albicans (Enjalbert et al., 2006; Enjalbert et al., 2003). There was only minor heterogeneity in the cellular response to CSP (Figure 3B).

Higher concentrations of CSP and longer time courses might induce a bifurcation similar to FCZ, leading to the emergence of Sd cells characterized by high expression of heat-shock protein coding genes. Since HSP70 was identified as a strong marker of the Sd response, we analyzed microscopy images of our GFP-tagged HSP70 cells exposed to CSP or FCZ for 2 and 3 days. Although a significant fraction of FCZ-treated cells highly expressed HSP70 at days 2 and 3, confirming the presence of Sd cells, no similar trend was observed in CSP-treated cells at either time point (Appendix 2—figure 1).

In conclusion, cells treated with CSP closely co-clustered with FCZ Rd response cells, but there was no evidence of CSP cells in the Sd response at day 3. The similarity in the early response is intriguing, considering that CSP and FCZ have distinct mechanisms of action (i.e., disrupting cell wall biosynthesis versus disrupting ergosterol biosynthesis). It will be interesting to determine the universality of the Rd response in future antifungal profiling efforts and the potential transition to the proliferative tolerant Sd state.

RAPA survivors exhibit the highest expression of hyperosmotic stress and lowest expression of RPs

The TOR (target of RAPA) pathway controls many processes including protein translation, autophagy, apoptosis, and cell growth in response to nutrient availability (Bastidas et al., 2009; Cruz et al., 2001; Heitman et al., 1991; Schmelzle and Hall, 2000). Although C. albicans is sensitive to RAPA (Baker et al., 1978), it has weak antifungal properties when used alone (Baker et al., 1978; Tong et al., 2021). However, RAPA treatment in combination with azoles, including FCZ, exhibits a synergistic effect (Tong et al., 2021), eliminating drug tolerance (Rosenberg et al., 2018). We sc-profiled C. albicans populations exposed to a concentration of RAPA that does not visibly impact growth (0.5 ng/ml) but is known to inhibit FCZ tolerance (Rosenberg et al., 2018).

Cells treated with RAPA primarily localize to cluster 2-pink and, to a lesser extent, cluster 4-turquoise in the Sd state (Figure 3A and B Appendix 2—figure 1). We previously noted that Sd FCZ cells exhibit high activation of heat-shock stress response genes. A similar trend is observed in RAPA-treated cells, but they mostly display a pronounced hyperosmotic stress response (Appendix 2—figure 1 ). In fact, this elevated hyperosmotic stress response is present in all RAPA cells, regardless of their cluster classification (Appendix 2—figure 1 , light blue curves). Specifically, RAPA cells express high levels of the sodium pump ENA2 and the glycerol-3-phosphate dehydrogenase GPD2 (Appendix 2—figure 1 ). Both genes are orthologs of hyperosmotic stress-responsive genes in S. cerevisiae and are induced in C. albicans in response to hyperosmotic stress (Enjalbert et al., 2003), indicating that hyperosmotic stress gene activation is unique to RAPA treatment.

Although RAPA-enriched regions initially seem to have moderate to low RP levels (Figure 3F and J), closer examination reveals that RAPA cells have the lowest RP expression compared to cells from any other condition (Appendix 2—figure 1). This is expected, as RAPA’s inhibition of TOR signaling is known to decrease RP gene expression (Powers and Walter, 1999). Consequently, no RAPA cells are found in the Rd state, which is characterized by high RP. When RAPA and FCZ are used together, drug tolerance is eliminated, but drug resistance remains unaffected. Our findings suggest that this might happen due to the inhibition of the Rd response, preventing cells from transitioning to the tolerant Sd state. Interestingly, RAPA specifically induces the expression of hyperosmotic stress-related genes. Past research has shown that the Tor1 pathway regulates phase separation and the biophysical properties of the cytoplasm by modulating macromolecular crowding (Delarue et al., 2018). The hyperosmotic stress response to RAPA exposure could be a result of disrupted cytoplasmic biophysical properties, leading to altered osmotic balance in cells.

Appendix 2—figure 1
Download asset Open asset
Single cell analysies of additional anti-fungal drugs.

(A, B, D–F) Uniform Manifold Approximation and Projection (UMAP) embedding from Figure 3A. Color bar depicts: (A) density of untreated (UT) cells, (B) expression score for the glycolysis signature. (C) Empirical cumulative distribution (ECDF) of the glycolysis signature scores for each cell in the different conditions (drug_day, color); Kolmogorov-Smirnoff test, p < 0.001. (D) Density of caspofungin (CSP) cells at day 2. (E) Histograms depicting expression of heat-shock protein 70 (HSP70) in (i) CSP or (ii) fluconazole (FCZ)-treated cells after 2 days (green) 29 and 3 days (blue) exposure. (F) Density of rapamycin (RAPA) cells after 2 days of exposure. Cells with high expression of HSP70 fall to the right of the dashed red line. There is a statistically significant depletion of cells highly expressing HSP70 (a marker of Sd; KM test, p < 0.01). (G) Score for the hyperosmotic stress signature mapped onto the UMAP embedding from Figure 3A. (H) ECDF of the hyperosmotic stress (hyperosmosticS) signature score for each cell in clusters 2 (top) and 4 (bottom) for the different conditions (drug_day, color). (I) The heatmap depicts genes (rows) in the hyperosmotic stress signature with significant variability and consistent expression across cells (columns). We used MAGIC to impute expression levels (z-score, color bar). Cells are linearly ordered by the overall magnitude of gene expression (mean gene rank). ‘Clusters’ (top) indicates the Leiden cluster of origin for each cell and ‘drug_day’ indicates the condition for each cell. (J) ECDF of signature scores for each cell under the different conditions (color) (ks-test, Kolmogorov-Smirnoff test).

Appendix 2—figure 2
Download asset Open asset
Differentially expressed genes and pathways associated with comets.

(A) Genes overexpressed in the collection of 19 comets compared to other cells in main clusters (Bayes factor > 2.5 and proportion of non-zero value > 0.2). Color bar indicates the mean expression across the cells in each cluster. The size of the dot indicates the proportion of cells which express that gene within each cluster. (B) Similar plot than (A) for genes differentially expressed in each comet cluster compared to others (Bayes factor > 3 and proportion of non-zero value > 0.2). (C) Gene ontology (GO) enrichment for genes identified as differentially expressed in comets compared to other clusters using pseudo-bulk differential expression analysis (60 genes; FDR < 0.1 ; Appendix 2—figure 2—source data 1, 2).The size of the dot is proportional to the number of genes in the list which overlap with the corresponding GO term. (D) Percent of reads assigned to genes in chromosome 2 for each cell classified in the 24 distinct clusters.

Appendix 2—figure 2—source data 1

Differentially expressed genes in cells classified in the comet clusters (pseudo-bulk samples) compared to other fluconazole (FCZ)-treated cells (pseudo-bulk samples).

baseMean = the average of the normalized counts taken over all pseudo-bulk samples; log2FoldChange = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange estimate; stat = Wald statistic; pvalue = Wald test p-value; padj = Benjamini-Hochberg adjusted p-value.

https://cdn.elifesciences.org/articles/81406/elife-81406-app2-fig2-data1-v2.xlsx
Appendix 2—figure 2—source data 2

Gene ontology terms enriched in genes differentially expressed between comets (pseudo-bulk) and other fluconazole (FCZ)-treated cells (pseudo-bulk) listed in Figure 3—table supplement 1A.

https://cdn.elifesciences.org/articles/81406/elife-81406-app2-fig2-data2-v2.xlsx

Data availability

Python/R code and data required for reproducibility is available through the Open Science Foundation (OSF) repository https://osf.io/5tpk3/ and associated github repository https://github.com/vdumeaux/sc-candida_paper (copy archived at Dumeaux, 2023). The raw and processed single-cell transcriptome and bulk RNA-seq is also available through NCBI's Gene Expression Omnibus with accession number GSE204903.

The following data sets were generated
    1. Dumeaux V
    2. Massahi S
    3. Bettauer V
    4. Khurdia S
    5. Omran RP
    6. Simpson S
    7. Xie JL
    8. Whiteway M
    9. Berman J
    10. Hallett MT
    (2022) NCBI Gene Expression Omnibus
    ID GSE204903. Candida albicans exhibits heterogeneous and adaptive cytoprotective responses to anti-fungal compounds.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE204903
    1. Dumeaux V
    (2022) Open Science Framework
    ID 5tpk3. sc-candida_paper.
    https://osf.io/5tpk3/

References

  95. Software
    1. R Development Core Team
    (2021) R: A language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
    http://www.r-project.org
  115. Book
    1. van Rossum G
    2. Drake FL
    (2009)
    Python 3 Reference Manual
    Scotts Valley, CA: CreateSpace.
    1. Wertheimer NB
    2. Stone N
    3. Berman J
    (2016) Ploidy dynamics and evolvability in fungi
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 371:20150461.
    https://doi.org/10.1098/rstb.2015.0461

Decision letter

  1. Antonis Rokas
    Reviewing Editor; Vanderbilt University, United States
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel
  3. David Shore
    Reviewer; University of Geneva, Switzerland

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Candida albicans exhibits heterogeneous and adaptive cytoprotective responses to anti-fungal compounds" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: David Shore (Reviewer #1).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

As you will see from the reviewers' individual reports, there was consensus that the manuscript and work reported are interesting but also had concerns about presentation and some of the controls used in the study. Thus, the manuscript will require revision, including potentially additional experiments, and one more round of review.

1) There is concern that there is significant overlap (and lack of clarity) with respect to data previously published (Battauer et al. 2020) and data that are new to this study. It will be imperative to clearly explain where the data in this new manuscript are coming from.

2) There was consensus among reviewers for a revised Introduction section that better sets up the rationale for this study (e.g., Reviewer #2, point #2 and Reviewer #3, point #1).

3) There was concern about the use of different time-points associated with controls vs. treatments (e.g., Reviewer #2, point #4) – addressing this point might necessitate collection of additional data.

Reviewer #1 (Recommendations for the authors):

Since the underlying biology here will be of most interest to researchers who are unlikely to be versed in the experimental and analytical methods used throughout, I think that the authors should make a better effort to help the reader begin to understand the bases of their analytical methods and their potential pitfalls. For example, they don't even define the ubiquitously applied "UMAP" method, much less explain in general how it works, what it reveals, and what its limits might be (though it is defined in the McInnes et al. reference title). I suspect that most readers will be unfamiliar with the details of single-cell transcriptomic analysis, which is complex and relatively recent. Μore science-related comments are as follows:

1. The significance of the dozen or so small "comet" clusters is unclear. On line 525 they are referred to as possible "trimeras" (What does that mean? Should it read "trisomes") "…with unstable genomes previously identified after FCZ treatment…which might have arisen as a response to amino-acid starvation." This latter point might be testable by increasing amino acid concentrations in the medium, which might then abolish the formation of some clusters.

2. In the Discussion the authors describe "bulk DNA-level profiling" of untreated cells and four FCZ treated populations (days 2, 3, 6, and 12). They state that attempts to assign the aberrations they observe to either the α or β response were inconclusive, but it's unclear how this would have even been possible to do. Please clarify. It is also stated in the text that variation levels peak at day 3, but it would seem to me to be day 2. Furthermore, on day 2 there would appear to be a nearly uniform increase in signal for all chromosomes shown, suggesting that there might be widespread trisomy. These points all need to be addressed and clarified.

3. The significance of the so-called "trajectory" in Figure 4A (black line) is unclear to me. From the text, it would appear to me that clusters 3 and 4 represent β-like cells on days 6 and 3, respectively. But what are the expression differences between these two groups? A similar question arises for clusters 1 and 2, which are said to be α-like cells. Does the black "trajectory" line imply a strict mode of transition between these 4 states? Why do α-like cells seem to be predominant when longer times are stated elsewhere to favor the β-like state?

4. The authors should state clearly how fast the cultures are growing (bulk measurement) at the time points they examine, and whether they have any way of knowing the growth rate as a function of transcriptome state and time of drug exposure. Related to this, do the authors imagine that α-like cells predominate early simply because they derive from those cells that had the highest levels of RP and ribosome biogenesis gene expression prior to treatment and thus might have been faster growing and/or more capable of escaping the early effects of drug treatment?

5. What is the actual evidence that β cells are derived from α cells through induction of RASTR in the latter population of cells? Related to his, could it be shown that α cells display a higher level of protein aggregation than β cells, perhaps by looking at RP fusions to GFP, versus HSP70-GFP?

6. Can the authors compare in more detail the transcription states of their β cells to the profile of RASTR cells (S. cerevisiae) described by Albert et al.? For example, are the targets of Hsf1 known in detail in Candida albicans, and if so, how well do they match with the up-regulated genes observed in β cells? Does Hsf1 also target genes encoding components of the proteasome (as in S. cerevisiae)?

7. Furthermore, can the authors show that Ifh1 is "condensed" in β cells, as seen clearly using Ifh1-GFP fusions in S. cerevisiae? In addition, what is known about Ifh1 targets in C. albicans (are they mostly RP genes?), and how well does this group overlap with the profiles of down-regulated genes in β cells? Perhaps the datasets are not robust enough to give this sort of information, but the issues should at least be addressed in the text by the authors.

8. The idea that a kind of persistent RASTR response promotes tolerance is very interesting. Perhaps the authors could test this idea further by asking whether inducing RASTR prior to drug treatment (with diazaborin or an RNA Pol I inhibitor) might strongly increase the fraction of tolerant cells.

9. Line 591: what is meant by "reinitiate translational machinery" with respect to the α cells? And how is this related to their persistence over time?

Reviewer #2 (Recommendations for the authors):

1. This manuscript emphasizes the lack of single-cell transcriptomics in C. albicans (and fungi broadly), although there does appear to be other work published in this area (Dohn et al. 2021 – which also includes antifungal treatment). Relationship to other work in this area should be more clearly addressed, and differences in findings with regards to antifungal treatment should also be addressed.

2. The introduction to this manuscript is framed around antifungal drug tolerance, and throughout this research, antifungal tolerance is highlighted as a central research question. However, it is not clear how the experimental design of this work addresses antifungal tolerance in any way. Cells are treated with antifungal drugs and subjected to transcriptomic analysis, but this does not represent a distinct 'tolerant' population of fungi that are being analyzed. This either needs to be much more clearly explained and justified, or more likely, the framing of this work needs to be substantially re-assessed.

3. The authors make reference to another study of theirs (Battauer et al. 2020) and suggest that the transcriptomic profiles reported in this work combine profiles from their previous work. This needs to be much more clearly explained. Have parts of this work already been previously published? How is this analysis unique and novel? The title of the previous manuscript seems very similar to that of this manuscript suggesting there might be a substantial overlap between these works.

4. The major analysis in this work compares untreated cells in log phase, to antifungal drug treated cells grown for 2-3 days in antifungals. It is not clear why untreated cells were not grown for the same duration of time as drug-treated cells, which would certainly alter the findings and analysis. It is thus unclear if the transcriptomic responses described in this manuscript truly represent the consequence of drug treatment itself, or are also influenced by the growth state of the cells (which are in stationary phase after 2-3 days of drug treatment). A comparison to untreated cells in stationary phase would be a more appropriate comparison.

5. The work describes the identification of 184 transcripts on average in each cell, which seems like an incredibly small number. Is this in line with other similar single-cell transcriptomic analysis? Does this number of transcripts enable robust conclusions to be drawn on the transcriptional profile of individual cells?

6. While the paper is generally well written, it is written in an extremely technical manner with much methodological detail (as well as discussion) incorporated throughout the main Results section. Major Results sections also lack clear and concise conclusions to help readers understand and interpret the major findings. This diminishes the clarity of the work and makes the manuscript at times quite dense and difficult to fully interpret.

7. It is unclear why stress response pathways are being assessed in untreated cells without any stress exposure. Should this be analyzed in the drug-treated cells instead?

8. In the growth curves in Figure 3 Supp 1a, it seems that both untreated and drug-treated cells take 2-3 days to reach stationary phase. This seems very unusual for cells that replicate quite rapidly under many growth conditions. Is this due to a nutrient-limited media or some other explanation?

9. It seems surprising that extremely different antifungals (cidal vs. static, different cell targets) elicit a very similar transcriptional response based on this analysis (line 345-348). Can this be explained?

10. Section 8-9 on expanding the analysis to day 6 was difficult to interpret in terms of the rationale for the experimental design and how the findings can be interpreted. It is also unclear how solid agar media-based tolerance assays can be extrapolated to liquid media growth assays with drugs, as these are very different conditions. It is also unclear if there is a day 6 untreated control that is being assessed.

a. The microscopy in associated Figure 4 K is very difficult to see and lacks scale bars.

Reviewer #3 (Recommendations for the authors):

Below I provide a few suggestions for how the authors may improve the manuscript, but overall I found it well done.

1. First, it is my understanding that tolerance is defined as survival, but not growth, after exposure above the MIC for drugs that normally would kill most of the population of cells (for non-static drugs). This is the definition put forward in Balaban et al. ("Definitions and guidelines for research on antibiotic persistence"; Nat Rev Microbiol; 2019). I would have called what the authors describe as phenotypic resistance, i.e., a sub-population of cells that can continue to grow after exposure above the MIC, and a non-genetically heritable state (i.e., readily reversible). This is a known phenomenon for some antimicrobials/species combinations, as the authors describe. I think it would be good to standardize usage of the word tolerance through the manuscript, or the authors can better explain how the phenomenon they observe is consistent with Balaban et al.'s definition of tolerance. This is an important semantics issue. Second, the drug concentrations used were at or well-below their MICs. This raises the question of whether the authors are studying the clinical phenomenon that they describe as tolerance (i.e., above the MIC) when the experiment was conducted below the MIC. While I understand their reasoning for using low drug concentrations, and I think there is still something to be learned at these concentrations, I think the authors should better contextualize the clinical relevance of their findings with the caveat that the experimental results were probably collected at sub-clinical concentrations. Admittedly, drug concentration is dynamic during treatment, so at some point the concentration in vivo likely passes through the author's choice of in vitro concentrations. It would be interesting to know how much their results are invariant to the drug concentration chosen.

2. The author's interpretation of the "comets" is questionable. On L244 they say, "This pattern suggests that the small set of cells from each comet have strong transcriptional similarity but each such comet is transcriptionally distinct from the other comets." I do not see how the cells within comets share "strong" transcriptional similarity because they are still spread out. Rather, the comet tails could be interpreted as lineages of descent, wherein each cell becomes more distinct (along some transcriptional axis) than its ancestor. This pattern is commonly observed in genetic data plotted with PCA, which is somewhat related to what is being shown here. Therefore, an alternative interpretation is that some lineages are moving toward a transcriptionally distinct state and leaving behind progeny along their trajectory, especially given that the experimental setup may not quickly remove ancestral cells in the two day exposure. Another explanation is that the comets are composed of a different (minority) morphology, such as filamentous growth. Another (less?) plausible interpretation is that these are noise clusters due to the inherent stochasticity involved in single-cell analyses. I see no reason why these explanations aren't equally reasonable to the author's explanation. Perhaps the authors can incorporate these alternative into their interpretations, but (at a minimum) there needs to be more justification given for the claim of "strong" intra-cluster transcriptional similarity.

3. L434: I am skeptical of the author's interpretation here about the relative fitness of α and β states. It could also suggest switching rates between phenotypic states are different between α and β populations. Also in this paragraph (L439), the authors use OD600 to show a higher growth rate for β over untreated. However, β also reaches a higher yield, which is unexplained. Having done many growth rate experiments with OD600 data, I am always skeptical of drawing large inferences about fitness from OD600. A better test of fitness is a competition experiment. I suggest the authors remove or downplay their conclusions about relative fitness here. Also, why isn't OD600 plotted in log-space (log(OD600)) if the goal is to see the difference in growth rate?

4. The authors did not seem to analyze feature importance in Leiden clustering, preferring to look at cluster associations (z-scores) instead. Leiden clustering is not guaranteed to be optimal, so some understanding of cluster stability would have been useful. That is, if inputs changed somewhat (e.g., subsampling or resampling), would the identified clusters have been similar? This is particularly important given that the UMAP1 versus UMAP2 clusters look like blobs that are split along arbitrary axes. It seems doubtful these clusters are stable, yet the whole manuscript analyzes them as though the blobs are discretized correctly.

5. I have some concerns about the comparison between untreated and treated at the same (or similar) time points. If the transcript profile changes over time, as presumably it would, then it is reasonable to expect that the treated and untreated cells are at different equivalent times because their growth rates are different. I have a similar concern with comparing different time points (days 3 versus 6) because fresh medium (YPD) was used in this experiment. How much are transcript abundances static over time? This should minimally be discussed further in the manuscript. Also, the conclusions in L472-475 and L579-592 may need this caveat mentioned.

https://doi.org/10.7554/eLife.81406.sa1

Author response

Essential revisions:

1) There is concern that there is significant overlap (and lack of clarity) with respect to data previously published (Battauer et al. 2020) and data that are new to this study. It will be imperative to clearly explain where the data in this new manuscript are coming from.

We understand this concern and we apologize for the lack of clarity in our previous version of this manuscript. We have adjusted the text in the introduction (and early Results section) to better explain this relationship. Briefly, we originally released the Bettauer et al. BioRxiv in 2020 (https://doi.org/10.1101/2020.01.21.914549), although we were unsatisfied with the descriptive nature of the single cell results, which primarily focused on the identification of the α and β responses. Our rationale was that two senior graduate students were set to finish their degrees, and we thought that the underlying technical contribution alone towards fungal single cell approaches was in and of itself sufficient to form a preprint and useful for the community (several other fungal groups had expressed interest in using our technology and we did not have the capacity ourselves to assist them directly).

However, we continued to do additional validation work (the fluorescence microscopy, CRISPR-mediated knockouts, DNA-sequencing and additional single cell RNA-sequencing) and this new data gave us deeper insight into the nature of the α and β responses. Bettauer et al. (2020) bioRxiv was not published elsewhere. Instead this data was merged with new data into this current effort.

Lastly, when we prepared for the first submission to eLife, we attempted to modify the Bettauer et al. (2020) preprint, but bioRxiv was tardy in re-opening the manuscript (we waited several weeks). We felt we could not delay submission any longer and submitted a new preprint instead: Dumeaux et al. (2022) bioRxiv (https://doi.org/10.1101/2022.07.20.500774). This data and analyses in this preprint are supersets of what appeared in Bettauer et al. (2020). We have asked that bioRxiv link this 2022 preprint to the original 2020 preprint. [Changes located at: L102-109; L116; L146-148].

2) There was consensus among reviewers for a revised Introduction section that better sets up the rationale for this study (e.g., Reviewer #2, point #2 and Reviewer #3, point #1).

We have re-worked the entire Introduction addressing specific comments from the Reviewers. We have provided more background and clearer definitions regarding fungal drug tolerance [Changes L53-61]. We have also sought to better explain the rationale for aspects of our study throughout this Response to Reviewers and in the manuscript (for example, for UT cells L184-187). We also moved a comparison of UT cells to the FCZ day 2 and 3 cells to Appendix 2 [Changes L 822-833]; this short comparison is not central to the main findings of the manuscript and has caused confusion as to our primary goal. We really only performed this comparison as a sort of “sanity” check to see if the single cell (sc)-profiling would detect differences between UT and drug treated cells.

3) There was concern about the use of different time-points associated with controls vs. treatments (e.g., Reviewer #2, point #4) – addressing this point might necessitate collection of additional data.

We fully agree with the concerns raised by the reviewers, and we have sought to better clarify throughout the manuscript our choice for controls. Our presentation within the previous version of the manuscript led to some confusion and we realize that we did not fully explain nor justify our choices in these dimensions. We address these issues below in response to specific questions from the reviewers and new data which was generated for this re-submission. [Changes L142-149 and several other sentences throughout the manuscript.]

A general note regarding this re-submission: we realized that the use of α and β to describe the two central subpopulations might be confusing to the C. albicans community, since α and a are used to describe the two mating types in this organism e.g. a/α, α/α, a/a. Globally now we refer to the α subtype as Rd (Ribo-dominant) and the β subtype as Sd (Stress-dominant) to avoid causing confusion in the literature.

Reviewer #1 (Recommendations for the authors):

Since the underlying biology here will be of most interest to researchers who are unlikely to be versed in the experimental and analytical methods used throughout, I think that the authors should make a better effort to help the reader begin to understand the bases of their analytical methods and their potential pitfalls. For example, they don't even define the ubiquitously applied "UMAP" method, much less explain in general how it works, what it reveals, and what its limits might be (though it is defined in the McInnes et al. reference title). I suspect that most readers will be unfamiliar with the details of single-cell transcriptomic analysis, which is complex and relatively recent.

We have revisited all aspects of the manuscript and inserted intuitive descriptions of the different methods used here. This includes the concepts of scVI UMAP embeddings, Leiden clusterings, pathway enrichment scores and trajectory analysis [L250-257, L368-373, L699-707]. We detail the location of additional changes below.

Μore science-related comments are as follows:

1. The significance of the dozen or so small "comet" clusters is unclear. On line 525 they are referred to as possible "trimeras" (What does that mean? Should it read "trisomes") "…with unstable genomes previously identified after FCZ treatment…which might have arisen as a response to amino-acid starvation." This latter point might be testable by increasing amino acid concentrations in the medium, which might then abolish the formation of some clusters.

We begin by noting that this section on “comet” clusters has been moved to Appendix 2 Supplemental Results as we feel that these observations are not critical to the manuscript and remain largely interesting observations [Changes L893-956]. The section has been completely re-written and clarification of the trimeras term is provided. Specifically, we were referring to our previous study in Harrison et al. where we observed three connected cells composed of a mother, daughter and granddaughter bud, termed trimeras, were observed within hours of FCZ exposure and were shown to give rise to aneuploidies at high frequency. This concept is not strictly necessary for this manuscript except to say that we do observe evidence of the associated aneuploidies [Changes L837-897].

It is an interesting idea to test if increased concentrations in the medium would ablate the appearance of such comets, and we certainly plan to explore this in the future. It is, however, somewhat tangential to the current effort and due to cost issues we haven’t explored this further to date. Disappearance of the comets might indicate that the GCN4-mediated amino acid starvation pathway is central to the formation of aneuploidies and/or provide insight into the rapid evolution of resistant micro-colonies post FCZ exposure. This in itself would constitute a significant manuscript.

Please also see our response to Reviewer #3, Comment 3 which touches on similar issues.

2. In the Discussion the authors describe "bulk DNA-level profiling" of untreated cells and four FCZ treated populations (days 2, 3, 6, and 12). They state that attempts to assign the aberrations they observe to either the α or β response were inconclusive, but it's unclear how this would have even been possible to do. Please clarify.

Before proceeding, we mention here that we decided in light of your comment to remove this paragraph from the Discussion section. Originally we intended to indicate that we did look at issues of chromosomal aberrations, which are common to C. albicans especially under stress, but our bulk DNA profiling was inconclusive. The sum of our data, analysis and validation experiments really require a full independent manuscript and the brief mention we make here is likely to only cause confusion.

To answer your technical question: The genes within each observed amplified genomic region were identified. We then examined whether those genes were simultaneously expressed more in either the α/Rd or β/Sd subpopulations according to the single-cell transcriptomics. The hope here was to find that a set of amplified genomic regions which were statistically enriched for genes that were over-expressed in exactly one of the α or β subpopulations. Since both the α/Rd and β/Sd subpopulations could “contribute” amplified genomic regions (and therefore over-expressed transcripts for their corresponding genes), we could perhaps assign some genomic regions uniquely to either the α or β subpopulations. We weren’t however able to conclude anything with statistical rigor here.

It is also stated in the text that variation levels peak at day 3, but it would seem to me to be day 2.

Yes indeed. This was a typo.

Furthermore, on day 2 there would appear to be a nearly uniform increase in signal for all chromosomes shown, suggesting that there might be widespread trisomy. These points all need to be addressed and clarified.

We agree. Note however the trimera occurrences described by Harrison et al. are not the only mechanism by which instability is observed in C. albicans in response to (drug) stress. In particular, the seminal papers by Ford et al. (doi: 10.7554/eLife.00662), Selmecki et al. (2009) and Todd and Selmecki (2020) (cited in the Introduction) describe other mechanisms for aneuploidy and loss across C. albicans. This includes recombination/amplification events involving specific-sequence motifs found across the C. albicans chromosomes. What we observe here may be due to other mechanisms besides the trisomy. For example, we did perform some ad hoc analysis comparing the predicted sites from Todd/Selmecki with regions of the C. albicans genome that were amplified in this study. Straightforward statistical approaches (based on the hypergeometric tests for enrichment) suggest that there is significant (albeit imperfect) overlap with their predictions. We also compared our regions of amplification with those of Ford et al. (see above) and again found significant, albeit imperfect, overlap. This suggests to us that there may be multiple mechanisms “controlling” genome instability in the C. albicans response to stress. The strength of amplification (black ticks) for chromosomes 4, 5, 6 and 7 are dramatic compared to other chromosomes (which appear closer to the pattern observed in the untreated samples), with only a few sporadic regions of amplification or loss. We have added to the text to clarify this analysis.

3. The significance of the so-called "trajectory" in Figure 4A (black line) is unclear to me. From the text, it would appear to me that clusters 3 and 4 represent β-like cells on days 6 and 3, respectively. But what are the expression differences between these two groups? A similar question arises for clusters 1 and 2, which are said to be α-like cells. Does the black "trajectory" line imply a strict mode of transition between these 4 states? Why do α-like cells seem to be predominant when longer times are stated elsewhere to favor the β-like state?

We note that we have re-written these paragraphs entirely to improve their readability. Changes Figure 4 is now Figure 5 pg 16. Changes in text L334-373. Discussion about the trajectory and meaning: [L368-373]

Indeed, clusters 3 and 4 represent β/Sd cells on days 6 and 3, respectively. First, we should emphasize that there are quite a few genes that are differentially expressed between β/Sd cells on days 6 and 3. The main differences are highlighted on L354-367 and Figure 3G.

In contrast, there are no significant robust difference for α/Rd day 3 versus day 6 (cluster stability analysis, Figure 5B).

Trajectory analysis based on single-cell RNA profiles allowed us to order cells along a “differentiation” trajectory. We now inserted intuitive descriptions for trajectory analysis in the manuscript; the essential idea is to consider the transcriptional profile of a cell as a sort of single still image from an animation. If the trajectory analysis (based on a tool called Slingshot) can order the still images to make a complete animation (denoted by the black lines), then this is considered one type of evidence that, when one population is adjacent to a second population (black edge), then the first population is derived from, or at least related to, the second population. In the end, the trajectory helps to confirm that β/Sd day 6 cells have evolved from β/Sd day 3 cells. For α/Rd, the lack of any robust transcriptional changes would already suggest that day 6 cells evolved from day 3 cells.

We think that there may be some confusion w.r.t. the last comment, as α/Rd cells are predominant at day 3 but β/Sd cells are predominant (or at least equal) at day 6. We speculate as to why β/Sd become more predominant on lines L364-367 and in the Discussion L393-399, L413-417.

4. The authors should state clearly how fast the cultures are growing (bulk measurement) at the time points they examine, and whether they have any way of knowing the growth rate as a function of transcriptome state and time of drug exposure. Related to this, do the authors imagine that α-like cells predominate early simply because they derive from those cells that had the highest levels of RP and ribosome biogenesis gene expression prior to treatment and thus might have been faster growing and/or more capable of escaping the early effects of drug treatment?

First, please note that we re-did the growth assays of cultures subjected to the different drugs at the relevant time points (Figures 3 figure supplemental 1B, Figure 5I). The cultures harvested for the untreated analysis were certainly in log phase; the growth curves suggest that the samples for all other time points are in stationary phase.

We do not find any evidence of a subpopulation of untreated cells that had higher RP and ribosome biogenesis expression and therefore we do not think that this is why the α/Rd cells predominate early. Our data suggests that the α/Rd cells only appear post treatment, suggestive of a “quick” acute response to drug exposure causing the dysfunction of the translational machinery.

5. What is the actual evidence that β cells are derived from α cells through induction of RASTR in the latter population of cells? Related to his, could it be shown that α cells display a higher level of protein aggregation than β cells, perhaps by looking at RP fusions to GFP, versus HSP70-GFP?

6. Can the authors compare in more detail the transcription states of their β cells to the profile of RASTR cells (S. cerevisiae) described by Albert et al.? For example, are the targets of Hsf1 known in detail in Candida albicans, and if so, how well do they match with the up-regulated genes observed in β cells? Does Hsf1 also target genes encoding components of the proteasome (as in S. cerevisiae)?

Our primary sources of evidence are the transcriptional profiles of α/Rd and β/Sd now presented in more detail in Section 5 of the Results. We expanded our analyses and have re-written some text to better express our core findings [changes L305-329]. One source of evidence for this relationship originates from analyses of the list of RASTR genes (and its associated processes, Table 1) from the literature and how they are expressed in either the α/Rd and β/Sd subpopulations. To investigate your questions in more detail, we identified C. albicans orthologs of S. cerevisiae RASTR genes (gene list supplied by B. Albert) together with a list of constitutive Hsf1 target genes identified in C. albicans (Leach et al., 2016). These genes, and their direction of expression, are now listed in Table 1A and B. Note for example that β/Sd cells exhibit decreased gene expression of ribosome processing (RP) genes (Figure 3E-F), consistent with the RASTR list in Table 1A (from B. Albert).

We also investigated whether β/Sd cells have enhanced expression of genes that are upregulated in RASTR by developing a signature score presented in Figure 4A (Albert list). Figure 4A shows that the right side of the UMAP, which is the location of the β/Sd cells, more highly express the genes listed as upregulated by Albert at al. (protein folding, proteolysis, glucose and pyruvate metabolic process, and others). Similarly, Figure 4B shows the expression of Hsf1 target genes from Leach et al. Again, the β/Sd cells show heightened expression of the constitutive targets of Hsf1. This comparison alone is not sufficient to argue that β/Sd cells are derived from the α/Rd cells. We can only hypothesize that β/Sd cells, which express all the hallmarks of RASTR activation, would survive and assume a transcriptional profile consistent with RASTR. At this time point, we do not know how we could track individual cells (e.g. via barcodes) and perform single-cell profiling to establish this dynamic behavior. The reviewer is correct to question our implication of causality. We have modified the text in the last section of the Results and Discussion to clarify this.

The observation that α/Rd cells predominate early (day 2) but there is an increase in the relative fraction of β/Sd cells at day 3 and day 6 also at least suggests a model whereby RASTR is “saving” α/Rd cells and sending them to β/Sd heaven.

At this stage, we do not know how to directly test for the relative degree of protein aggregation of α/Rd versus β/Sd cells. We had difficulties to establish a good dual fluorescent strain to measure both aggregation in the Rd and Sd cells simultaneously. We do observe evidence of aggregation, but we are unable to strictly limit this to α/Rd and not β/Sd cells. We do want to stress here that α/Rd transcriptional profiles are very distinct with overexpression of many genes almost exclusively related to ribosome biogenesis, rRNA processing and protein translation [L285-290; Figure 3 —figure supplement 2B; Figure 3 – source data 1,2].

7. Furthermore, can the authors show that Ifh1 is "condensed" in β cells, as seen clearly using Ifh1-GFP fusions in S. cerevisiae? In addition, what is known about Ifh1 targets in C. albicans (are they mostly RP genes?), and how well does this group overlap with the profiles of down-regulated genes in β cells? Perhaps the datasets are not robust enough to give this sort of information, but the issues should at least be addressed in the text by the authors.

These are excellent suggestions. Unfortunately, the Ifh1 gene is not well represented in our sc-transcriptomics data, and the meager detected expression is not significantly different between the α/Rd and β/Sd subpopulation. We focused instead on the expression profile of targets of Ifh1.

– Wade et al. (PMID: 15616568) state only that the vast majority of targets are RP genes. As per Figure 3F, we know that the RP genes are high in the α/Rd subpopulation and very low in the β/Sd subpopulation, consistent with the RASTR sensing mechanism.

– We examined Ifh1 targets provided by PathoYeastract (C.albicans http://yeastract-plus.org/pathoyeastract/calbicans/index.php) (n = 148 genes incl. 41 putative genes and 64 RP-encoding genes). Among those, 103 are detected in our profiles. We selected the most variable genes across cells (n=74) and plotted their expression in the heatmap in Author response image 1. Most targets are more highly expressed in α/Rd cells (darkpink cluster1) than in β/Sd cells (turquoise cluster4): this includes many RP-coding genes as well as YST1 (ribosome associated protein gene), ARF2 (stress response) and UBI3 (alias RPS31). The rest of the IFH1P targets (n = 10 at the bottom of heatmap) are most highly expressed in cluster (green cluster 3) enriched for untreated cells; they are not differentially expressed between the α/Rd and β/Sd subpopulations except for CCP1 (stress response), ERG251 (ergosterol biosynthesis), and ADE17 (adenine biosynthesis). These genes are present in the list of genes differentially expressed between α/Rd and β/Sd provided in Figure 3 – source data 1.

Overall, this analysis further confirms that RP-encoding genes are highly expressed in α/Rd cells (and lowly expressed in β/Sd cells). Since this is well established in the manuscript already, we choose to restrict this figure to Author response image 1 only. We have added text to discuss Ifh1. [Changes L326-329]

Author response image 1
Download asset Open asset
The heatmap depicts the expression of the most variable IFH1 targets (rows) across cells (columns).

We used MAGIC to impute expression levels (z-score, color bar). Cells are linearly ordered by the overall magnitude of gene expression (mean gene rank). “Clusters” (top) indicate in which cluster the cell is assigned. Recall that cluster 1 (red) corresponds to α/Rd and cluster 4 (turquoise) corresponds to β/Sd. Cluster 3 corresponds to untreated cells.

8. The idea that a kind of persistent RASTR response promotes tolerance is very interesting. Perhaps the authors could test this idea further by asking whether inducing RASTR prior to drug treatment (with diazaborin or an RNA Pol I inhibitor) might strongly increase the fraction of tolerant cells.

This is an excellent idea and we very much intend to do this and a series of additional supporting experiments with appropriate controls. After extensive discussion with our team, we judged that a proper in depth investigation of this would require a very significant and expensive effort. Although the results would add more mechanistic understanding/validation to this manuscript, it is currently not financially feasible for my lab and would delay publication of this already lengthy manuscript by another year. We truly hope that the reviewer understands our constraints and the costs associated e.g. with the sc-profiling.

9. Line 591: what is meant by "reinitiate translational machinery" with respect to the α cells? And how is this related to their persistence over time?

Thank you for pointing this out – we have removed this idea from the manuscript.

Reviewer #2 (Recommendations for the authors):

1. This manuscript emphasizes the lack of single-cell transcriptomics in C. albicans (and fungi broadly), although there does appear to be other work published in this area (Dohn et al. 2021 – which also includes antifungal treatment). Relationship to other work in this area should be more clearly addressed, and differences in findings with regards to antifungal treatment should also be addressed.

We have adjusted the Introduction to better cover the results of Dohn et al. However, we would like to note that the original bioRxiv (2020) version of this paper appeared well before Dohn et al. in 2021. It is an unfortunate oversight that Dohn et al. does not cite our effort as microbial DROP-seq (mDROP-seq) has similarities to our published protocol. Nevertheless, there are many differences in the tools and parameters between our two efforts with respect to the quality control and pre-processing of sc-transcriptomics data, yet the number of genes/transcripts per cell were quite comparable.

The results in Dohn et al. support our findings, including those reported in our original 2020 bioRxiv preprint. To summarize their results:

– A “mixed sample” of S. cerevisiae and C. albicans was DROP-seq’d. They could differentiate between these two fungi. This is reminiscent of the first DROP-seq experiments from Macosko et al. (2015) Cell where mixed mouse and human samples were used to establish the technical efficacy of their system. A nice experiment but independent of any of our results here.

– They examine S. cerevisiae subjected to heat shock. This is not related to our manuscript and is similar in structure to several previous single cell profiling efforts that looked at heat shock in S. cerevisiae, including the series of papers from the Gasch and later McClean labs. A nice experiment but we do not examine S. cerevisiae here so we cannot compare.

– They contrast control C. albicans at 1.5 hours and 3 hours post exposure to FCZ. We examine C. albicans populations at 2, 3, and 6 days post treatment. Although these are very different time scales, there is still some overlap between our differentially expressed genes (e.g. WH11 which is discussed in Bettauer et al., 2020). Their analysis of the cell cycle relies on older bulk transcriptional profiles from Spellman et al. By contrast, we combined several sources of evidence including cell cycle related information obtained from more modern transcriptomics data. It is hard to “align” our data with their data for these reasons. At 1.5-3 hours post treatment, Dohn et al. appear to have captured an early cellular stress response. The trajectory analyses mostly point to differences in cell cycle. [Changes located at L102-108]

2. The introduction to this manuscript is framed around antifungal drug tolerance, and throughout this research, antifungal tolerance is highlighted as a central research question. However it is not clear how the experimental design of this work addresses antifungal tolerance in any way. Cells are treated with antifungal drugs and subjected to transcriptomic analysis, but this does not represent a distinct 'tolerant' population of fungi that are being analyzed. This either needs to be much more clearly explained and justified, or more likely, the framing of this work needs to be substantially re-assessed.

We thank the reviewer for pointing out aspects of the introduction that were not clear. We indeed designed the study to focus on antifungal tolerance, by using drug concentrations at and just above the MIC (MIC of C. albicans lab strain SC5314 is 0.25-0.5 µg/ml fluconazole) and by analyzing cells at 48 h of drug exposure (a time at which the differences between tolerant and non-tolerant cell growth is evident) (Rosenberg et al. 2018; Gerstein et al. 2016). We now address this by highlighting several aspects of tolerance: First, when tolerance appears, some cells are dividing much more than others. The dividing cells are those exhibiting tolerance to fluconazole, which is a fungistatic drug. We now highlight this, as it is a major distinction between the definition of fluconazole tolerance in fungal cells being treated with a fungistatic drug and the definition of antibacterial tolerance in bacterial cells being treated with a bactericidal drug. The question we are asking is whether there are major distinct subpopulations of cells exhibiting different responses to drug within a single population of isogenic cells. The answer to this is clearly yes, especially with respect to fluconazole, and the two major types of cells are those with an α/Rd and those with a α/Rd and β/Sd response pattern. [Changes L53-73].

3. The authors make reference to another study of theirs (Battauer et al. 2020) and suggest that the transcriptomic profiles reported in this work combine profiles from their previous work. This needs to be much more clearly explained. Have parts of this work already been previously published? How is this analysis unique and novel? The title of the previous manuscript seems very similar to that of this manuscript suggesting there might be a substantial overlap between these works.

Please see our response to these questions above (Essential Revisions #1 and the changes located at L96-99; L109;-114; L137-141). Briefly here, the first preprint from the Hallett lab was never published beyond bioRxiv. Here, this manuscript uses that original bioRxiv data in addition to new data including: more replicates of the sc-transcriptional profiles, as well as microscopy, CRISPR knockouts, DNA-sequencing and disk assays. This manuscript goes significantly further than the first bioRxiv paper, which primarily detailed the technical aspects of the sc-RNAseq assay and identified the existence of the α/Rd and β/Sd survivor subpopulations after FCZ treatment. Several additional authors were added to this manuscript to assist us with the interpretation of the data (for example, the Berman team, which provided antifungal drug tolerance expertise and performed some of these additional studies).

4. The major analysis in this work compares untreated cells in log phase, to antifungal drug treated cells grown for 2-3 days in antifungals. It is not clear why untreated cells were not grown for the same duration of time as drug-treated cells, which would certainly alter the findings and analysis. It is thus unclear if the transcriptomic responses described in this manuscript truly represent the consequence of drug treatment itself, or are also influenced by the growth state of the cells (which are in stationary phase after 2-3 days of drug treatment). A comparison to untreated cells in stationary phase would be a more appropriate comparison.

There appears to be a misunderstanding. The untreated cells were in log phase, which allowed us to define cell-cycle-specific gene expression patterns. The tolerant cells exposed to drug were still growing, but much more slowly, and they continued to grow between days 2 and 3 of the experiment (Figure 3 —figure supplemental 1B). Between days 3-6, the cells were returned to medium without drug so that recovery could be documented, although in this case, because they had three days to recover, the cells were well past log phase when they were analyzed. It is true that many factors beyond the drug may affect heterogeneity. We think that the interesting observation in these results is that, even after 3 days of recovery, the α/Rd and β/Sd split persisted. In addition, we note that the α/Rd and β/Sd cell clusters are defined solely within each culture and neither by contrasting cultures between time points nor by comparing UT versus drug treatments.

We have made more clear in the Methods and in the manuscript in several places e.g. [L196, L509], that state of the cells at time of profiling for each drug/day. We also re-perform the growth curves (Figure 3 figure supplemental 1B and Figure 5I) to provide the reader with a more quantitative estimate of the state of these cells at time of profiling.

5. The work describes the identification of 184 transcripts on average in each cell, which seems like an incredibly small number. Is this in line with other similar single-cell transcriptomic analysis? Does this number of transcripts enable robust conclusions to be drawn on the transcriptional profile of individual cells?

Indeed, expressed in this manner, 184 transcripts per cell on average appears to be a very small number. However, we strongly caution the reviewer from concluding too much based only on the first moment (mean) of a distribution in this context. First, note that there is a very long right tail representing cells with far more transcripts; these cells strongly influence the downstream analyses.

Specifically, every sc-transcriptomics method begins with a quality control pipeline that selects “high quality” cells. The definition of high quality is based on approximately 10 underlying parameters that are determined by different types of ad hoc statistical analyses and specific bioinformatics tools (e.g. EmptyDrops software that seeks to determine thresholds of detected RNA that represent successful co-capture of a cell in a droplet versus sequencing of ambient RNA ubiquitously present in the suspension). There are many ways to increase the mean transcripts per cell at the expense of throwing out some sparse cells; however these sparse cells may contribute useful information to the overall patterns in the data. For example, key markers of the α/Rd and β/Sd classes may be expressed in these “sparse cells” and assist to delineate between subpopulations.

This leads to the question: If you cannot maximize both the number of cells you sequence and the number of transcripts you sequence per cell, which would you choose? It turns out that it is better to increase the number of cells (Zhang, Ntranos, Tse (2020) Nature Communications), even though this may seem counterintuitive at first. A standard approach in sc studies is to first use the sparse sc-data to identify subpopulations. Once identified, we can then pool the data within each subpopulation (ignore barcodes) and do traditional differential analysis. This is exactly what we did in this manuscript: we used the sparse single-cell profiles to discover the α/Rd and β/Sd subpopulations. We then used traditional bulk analysis by simply binning cells as either α/Rd or β/Sd (or neither) and then ignored the cell barcodes. This approach provided pseudo-bulk RNA-seq profiles of the α/Rd and β/Sd subpopulations that were quite robust. It also allowed us to drill down into the genes that are differentially expressed between the two subpopulations.

As we are the first to do sc sequencing in C. albicans at large scale, it is a challenge to establish comparisons with previous studies. A study in Saccharomyces performed 11 replicates of their samples with a commercial 10x Inc system (Jackson et al. (2020) eLife). Nevertheless, they had a very similar number of transcripts per cell to our data – we captured a maximum of 5956 transcripts/cell and they captured a maximum 5437 transcripts/cell. Note, however, that the lower bound of transcripts/cell is a function of many parameters and, if we repeated the DROP-seq nine more times, we could discard cells to raise our lower bound to an equivalent of 735 transcripts per cell reported by Jackson et al. At the time the experiments were done, each run on the Chromium 10x system was approximately 10x more expensive than a run using our system.

The numbers from Dohn et al. were almost identical to the numbers in this manuscript, although there are differences in many parameters: concentration of beads used, number of DROP-seq runs combined before sequencing, depth of sequencing, quality control parameters etc. Overall, they profile a smaller number of cells (however we consider multiple drugs and more time points). However, it does really come down to Zhang et al.’s observation above – it is more important to increase the number of cells analyzed, than to capture more transcripts per cell. Our goal was to show proof or principle for fungal sc-transcriptional profiling at a reasonable cost, and we think we have achieved that: we were able to identify and partially validate subpopulations that can be explored more deeply.

We also note that section 2 of the Results describes comparisons with bulk transcriptional profiling (Methods 5) and bulk differential expression analyses (Methods 9), which support the idea that the sc-profiling is largely detecting the same transcripts (~6700 in common, only 172 that differ).

In conclusion, whether 184 transcripts per cell (for our specific choice of thresholds and parameters) is or is not sufficient depends on the complexity of the samples, the types of questions that are investigated and the ability to independently validate the existence of subpopulations. Here, for FCZ treated cells for example, we only needed sc-profiles to be of sufficient quality to separate cells into the α/Rd and β/Sd subpopulations. Improvements to the technology would allow us to identify rarer subpopulations.

6. While the paper is generally well written, it is written in an extremely technical manner with much methodological detail (as well as discussion) incorporated throughout the main Results section. Major Results sections also lack clear and concise conclusions to help readers understand and interpret the major findings. This diminishes the clarity of the work and makes the manuscript at times quite dense and difficult to fully interpret.

We understand and have made an effort to ensure that each section contains intuitive, non-expert explanations of the main findings. We also moved some of the more technical and less relevant material to the appendices (e.g. discussion of comets, results for the CSP- and RAPA- treated cells) and reformulated some sections of the paper to clearly highlight when the text refers to hypotheses concerning the biological interpretation of our results. Our goal was to focus the manuscript primarily on the discussion of the Rd (α) and Sd (β) subpopulations.

Our response to your previous question on page 10 provides the line numbers in the new manuscript where these changes occur.

7. It is unclear why stress response pathways are being assessed in untreated cells without any stress exposure. Should this be analyzed in the drug-treated cells instead?

We have modified the text to only discuss pathways relevant to untreated (UT) cells in section 3 (which focuses on UT cells). In Figure 2 – source data 1 (list of gene signatures), we added a column entitled “UT relevant”, which indicates which of the 43 curated gene signatures is relevant for analyses associated with profiles of UT cells. Note however that there is evidence in the literature that some stress pathways do have (at least slight) differential expression across different cell cycle states. This includes, for example, heat shock proteins. This is reflected in the data where the signatures corresponding to these different stress pathways are also differentially expressed (for example, the oxidative stress response). [Changes L188-197]

8. In the growth curves in Figure 3 Supp 1a, it seems that both untreated and drug-treated cells take 2-3 days to reach stationary phase. This seems very unusual for cells that replicate quite rapidly under many growth conditions. Is this due to a nutrient-limited media or some other explanation?

We apologize for this confusion. Indeed, the cells were grown in YPD media so there should not have been any nutrient limitation per se. Although the underlying growth experiments had been repeated several times, we were concerned that there may have been some error in the measurement. Thus, we repeated these experiments and found different results. We now replaced these figures with the new data (Figure 3 figure supplemental 1B; Figure 5I); indeed the stationary phase is reached much earlier, and drug-treated cells exhibit slower growth as expected. We suspect that the confusion was due to a failure to account for the dynamic range of the spectrophotometer.

9. It seems surprising that extremely different antifungals (cidal vs. static, different cell targets) elicit a very similar transcriptional response based on this analysis (line 345-348). Can this be explained?

To be clear here, we observed that the majority of caspofungin (CSP) survivors have similar transcriptional profiles to FCZ cells in the α/Rd state; we did not see any CSP survivors with profiles resembling the β/Sd state, so there are some similarities and some differences between CSP, which is fungicidal in C. albicans and FCZ, which is fungistatic. Of note, because CSP is fungicidal, we used subinhibitory CSP concentrations (to ensure that a sufficient number of cells survived for profiling) and the CSP sc-profiling was only done at day 2. Nonetheless, we did examine cultures at later time points and with higher concentrations of CSP using microscopy. We always observed high levels of cell-cell aggregation which makes them inadequate for sc-profiling, yet did not observe any indication that CSP survivors of type β/Sd under any conditions (see eg Appendix 2 – figure 1E depicting Hsp70 fluorescence, marker of β/Sd cells, at day 2-3 for CSP-treated cells). In summary, perhaps cells at early time points post-treatment with low doses of CSP primarily cope with protein aggregation, a property shared with FCZ, but not with other tested compounds such as Rapamaycin.

We have added text to the Discussion section on the CSP response L428-435 and additional text regarding the RAPA response (which conversely only induces β/Sd cells L418-427).

10. Section 8-9 on expanding the analysis to day 6 was difficult to interpret in terms of the rationale for the experimental design and how the findings can be interpreted. It is also unclear how solid agar media-based tolerance assays can be extrapolated to liquid media growth assays with drugs, as these are very different conditions. It is also unclear if there is a day 6 untreated control that is being assessed.

Indeed, agar assays (disk diffusion assays, DDAs) and liquid assays (broth microdilution assays) do not give exactly the same results. This is because, in liquid cultures, cells grown in the presence of a drug can out-compete those cells that fail to grow in the drug. Nonetheless, when assayed at 48 hours, both disk assays and broth microdilution assays can report on tolerance (Rosenberg et al. 2018). Usually tolerance is measured at 48 hrs for cultures that grow at normal rates, and sometimes at 72 hrs for slow-going cultures. The reviewer is correct to question the meaning of the DDAs presented in Appendix 1 at day 6: whereas the cells harvested for sc-profiling at day 6 are return-to-growth without drug, the 6 day measurement in the DDA has not received replenished media. Nevertheless, the DDAS of Appendix 1 still indicates that the cultures do exhibit drug tolerance in the 2-3 day window relevant to our sc-profiling and that we do not observe drug resistance. We have added text in Appendix 1 to clarify this. [Changes L122-131 and in Appendix 1 L782-820].

a. The microscopy in associated Figure 4 K is very difficult to see and lacks scale bars.

We have decided to remove this panel; the two microscopy images were primarily there to show that we observed cells at the day 6 time point. The images do not advance the main findings in the manuscript significantly.

Reviewer #3 (Recommendations for the authors):

Below I provide a few suggestions for how the authors may improve the manuscript, but overall I found it well done.

1. First, it is my understanding that tolerance is defined as survival, but not growth, after exposure above the MIC for drugs that normally would kill most of the population of cells (for non-static drugs). This is the definition put forward in Balaban et al. ("Definitions and guidelines for research on antibiotic persistence"; Nat Rev Microbiol; 2019). I would have called what the authors describe as phenotypic resistance, i.e., a sub-population of cells that can continue to grow after exposure above the MIC, and a non-genetically heritable state (i.e., readily reversible). This is a known phenomenon for some antimicrobials/species combinations, as the authors describe. I think it would be good to standardize usage of the word tolerance through the manuscript, or the authors can better explain how the phenomenon they observe is consistent with Balaban et al.'s definition of tolerance. This is an important semantics issue.

Indeed we are discussing a phenotypic heterogeneity, (which is due to slow growth during constant exposure to a static drug); reviewers of prior work that characterized this phenomenon insisted that it be termed `tolerance’ because some older studies had referred to it in this way. Thus, although we suggested using a different term (perseverance) to distinguish this phenomenon from bacterial tolerance (which indeed is due to survival and not growth during periodic exposure to a cidal drug), editors of that paper (Rosenberg et al. 2018) agreed with the reviewers, and thus we are stuck with as suboptimal term which has different meanings in different contexts. For this manuscript, we tried to clarify this issue in the introduction (and to point out that the definition differs from that in bacteria treated periodically with cidal drugs). [Changes 62-73]

Second, the drug concentrations used were at or well-below their MICs. This raises the question of whether the authors are studying the clinical phenomenon that they describe as tolerance (i.e., above the MIC) when the experiment was conducted below the MIC. While I understand their reasoning for using low drug concentrations, and I think there is still something to be learned at these concentrations, I think the authors should better contextualize the clinical relevance of their findings with the caveat that the experimental results were probably collected at sub-clinical concentrations. Admittedly, drug concentration is dynamic during treatment, so at some point the concentration in vivo likely passes through the author's choice of in vitro concentrations. It would be interesting to know how much their results are invariant to the drug concentration chosen.

There appears to be some misunderstanding: the FCZ MIC50 of C. albicans SC5314 is 0.5-1 µg/ml. Here we used 1µg/ml fluconazole which is above the MIC (1-2x MIC50). We had previously stated ~1x MIC50 as this was an estimate based on MIC50 obtained from the literature. We have added a broth microdilution experiment to better estimate the FCZ MIC50 for the specific strain used in this study. The result now presented in Figure 3 —figure supplement 1A for FCZ at days 2 and 3, and Figure 5I for day 6.

Indeed, antifungal tolerance occurs at higher drug concentrations as well. However, antifungal tolerance is often very similar at drug concentrations ranging from just above the MIC to concentrations far above the MIC (e.g., Rosenberg et al. 2018, Nat. Comm.; Yang et al., mbio 2023; Todd et al., 2023). The reviewer is correct that the concentrations used here are sub-clinical and we agree that it would be very interesting to look at this question in a series of drug concentrations. This is a major challenge however in and of itself, and the cost of profiling/sequencing samples across a range of FCZ concentrations and time points with suitable controls is prohibitive, and would require substantial funding. In general, we agree it would be very interesting to better understand the dynamics of fungal drug tolerance in in vivo models relevant to the clinical setting and we hope to do so in the future.

2. The author's interpretation of the "comets" is questionable. On L244 they say, "This pattern suggests that the small set of cells from each comet have strong transcriptional similarity but each such comet is transcriptionally distinct from the other comets." I do not see how the cells within comets share "strong" transcriptional similarity because they are still spread out. Rather, the comet tails could be interpreted as lineages of descent, wherein each cell becomes more distinct (along some transcriptional axis) than its ancestor. This pattern is commonly observed in genetic data plotted with PCA, which is somewhat related to what is being shown here. Therefore, an alternative interpretation is that some lineages are moving toward a transcriptionally distinct state and leaving behind progeny along their trajectory, especially given that the experimental setup may not quickly remove ancestral cells in the two day exposure. Another explanation is that the comets are composed of a different (minority) morphology, such as filamentous growth. Another (less?) plausible interpretation is that these are noise clusters due to the inherent stochasticity involved in single-cell analyses. I see no reason why these explanations aren't equally reasonable to the author's explanation. Perhaps the authors can incorporate these alternative into their interpretations, but (at a minimum) there needs to be more justification given for the claim of "strong" intra-cluster transcriptional similarity.

We thank the reviewer for these insights. Because we do not consider the “comets” to be central to the main findings of the manuscript, and because we think that they potentially detract and distract from the central points, we decided to move the section on comets to the Appendix 2. In addition, we modified the text concerning comets to clarify the previous presentation [Changes L837-897]. Below we provide specific responses to the points raised above concerning interpretations of the comet data.

"This pattern suggests that the small set of cells from each comet have strong transcriptional similarity but each such comet is transcriptionally distinct from the other comets." I do not see how the cells within comets share "strong" transcriptional similarity because they are still spread out.

There are 19 smaller clusters or “comets” outside of the 5 major clusters analyzed in the main body of the manuscript. Each comet cluster is, in turn, comprised of a set of 7 to 126 cells. The cells within each comet are similar at the transcriptional level – the UMAP embedding procedure groups them together based only on their transcriptional profiles. However, the 19 different comets are not necessarily transcriptionally similar to one another, as noted, because “they are still spread out” and found in distinct clusters. [Changes 837-839]

Nevertheless, when we pooled all 19 comets and computed the list of all genes that were differentially expressed in the pooled comets relative to the 5 main clusters, we did find some interesting genes and pathways that suggest that these “outliers” are not simply technical artifacts of the DROP-seq profiling or the UMAP clustering algorithm. In particular, we discuss the ILV pathway which is over-represented in this group and its relationship to amino-acid starvation amongst other pathways revisited below.

We appreciate the original presentation of the analysis was dense and thus, we rewrote the section to better explain the use of the different tools and statistics. The comets are small and the sc-profiles are sparse, so analysis here really benefits from alternative methods. We have tried to better delineate these different sources of evidence [e.g. italic subtitles for each method on L844, 854, and 882]. Fortunately the different approaches largely highlight the same take home message, which we have tried to clarify in the presentation. [Changes L837-897]

We also removed the somewhat poetic language of “radiating from the centre of the UMAP”. We fully agree that this sort of “trajectory” is often reflective of some sort of model underlying the decomposition technique whether PCA (which is linear) or UMAP (which could be non-linear). This “linear radiation” first suggested to us that the comets might be some sort of technical artifact, but the differentially expressed genes and the pathway enrichment analysis suggests alternative explanations that have some support in the previous literature of C. albicans drug tolerance and its potential relationship with the genomic instability observed in C. albicans. [Changes 267-269]

We caution against interpreting the linear trajectories as true “evolutionary” trajectories. This is certainly one possibility, but the size of each comet and the current limitations of sc-transcriptional profiling would make a rigorous statistical assessment of this very difficult.

That the comets could be due to stochastic noise is another option, but we disagree that it is equally likely, given that both the differentially expressed gene analysis and the pathway analysis are tools underpinned by sound statistical models, and both provide evidence (at reasonable p-values or confidence intervals) that there is enrichment in specific functional processes. It is more difficult to estimate a formal confidence in the mapping of reads to chromosomal position (i.e., cluster 16 to chromosome 2) however it is at least visually convincing that ~38% of all reads from cells in cluster 16 mapped to this chromosome. We also note that comets are almost exclusively observed after FCZ treatment (and not after Rapamycin, for example). We have added a statement to the end of the section that all of our findings should be taken with caution due to the small sample sizes. Changes L889-890.

It is possible that some of these clusters correspond to filamentation, although (i) we filtered out large filamentous cells before profiling, (ii) filamentation signatures are not identified in the analysis, (iii) microscopy analysis of the cultures before filtering identified <1% of cells in a filamentous state. (Please see L148-149 and L492-493 in the manuscript.) It is possible that these cells could be representative of “small” germ tube cells, although we did not identify marker genes related to morphological switches.

This trifecta of analyses is why discussion regarding the comets remained in the manuscript – the enrichment of processes established to be central in fungal drug tolerance was far too high for this to be explained by “noise” (or the background models of any of these tools).

Our text was confusing (e.g. nature of intra- and inter-comet variability) and perhaps did not highlight why the extensive alternative analyses led us away from the idea that these were mere technical artifacts. However we agree that we are not able to fully characterize the comets at this point (e.g. exploring rigorously through DNA-seq what chromosomal aberrations exist and whether they line up well with the seminal work of Todd and Selmecki) so it remains descriptive at this point. We have minimized the discussion about DNA level aberrations (we also removed a paragraph about whole genome DNA profiling of the samples from the Discussion to simplify the manuscript).

3. L434: I am skeptical of the author's interpretation here about the relative fitness of α and β states. It could also suggest switching rates between phenotypic states are different between α and β populations. Also in this paragraph (L439), the authors use OD600 to show a higher growth rate for β over untreated. However, β also reaches a higher yield, which is unexplained. Having done many growth rate experiments with OD600 data, I am always skeptical of drawing large inferences about fitness from OD600. A better test of fitness is a competition experiment. I suggest the authors remove or downplay their conclusions about relative fitness here. Also, why isn't OD600 plotted in log-space (log(OD600)) if the goal is to see the difference in growth rate?

We fully agree with the skepticism regarding the relative fitness of α/Rd and β/Sd states, and these comments motivated us to revisit our group’s previous discussions as to what we can say regarding fitness and growth rates of the α/Rd and β populations. First, it is important to clarify that L439, which previously referenced Figure 4I (now Figure 5I) does not show β/Sd cells over untreated cells. The survivors at day 3 are a mixture of α/Rd and β/Sd cells. However, we did report that FCZ day 3 survivors grew faster after receiving new medium from days 3 to 6 compared to untreated cells, an observation that we struggled to explain. We decided to repeat the growth assays analysis of UT cells and FCZ day 3 survivors (Figure 5I) with care to ensure that both populations were diluted to the same OD600 initially (Methods 1C). (See response to Reviewer 2 Question 8). After correcting some prior issues with dilution and the dynamic range of the spectrometer, we now observe the more intuitive outcome where UT cells grow faster than FCZ survivors. Our previous manuscript stated that “At these concentrations, FCZ treatment slows growth relative to UT controls in the first days after exposure” which is now more visible with growth curves fitting what you would expect. We are sorry for this confusion.

Currently, we do not have a way to directly test the fitness of the α/Rd versus β/Sd survivors. We do comment however on the relative sizes of the α/Rd subpopulations compared to β/Sd subpopulations. In particular, from day 2 to day 6, we see that the β/Sd becomes more prevalent (L359-367). We do not have a mechanistic explanation but we offer at least two different models that are both consistent with our data as to why this may be the case (Changes L364-367 and L393-399, L413-417).

4. The authors did not seem to analyze feature importance in Leiden clustering, preferring to look at cluster associations (z-scores) instead. Leiden clustering is not guaranteed to be optimal, so some understanding of cluster stability would have been useful. That is, if inputs changed somewhat (e.g., subsampling or resampling), would the identified clusters have been similar? This is particularly important given that the UMAP1 versus UMAP2 clusters look like blobs that are split along arbitrary axes. It seems doubtful these clusters are stable, yet the whole manuscript analyzes them as though the blobs are discretized correctly.

This is an excellent point. We now include an updated analysis that incorporates measures of cluster stability throughout the manuscript. Before describing our changes and additional analyses, we want to note that the vast majority of single cell transcriptomic studies are performed with mammalian tissues. The diversity of cellular lineages and cell types (and extreme disruptions of normal physiological states in studies related to disease) induce UMAPs with distinct, non-overlapping clusters. We hypothesize here that this degree of separation will not be observed in a fungal population where the cells are isogenic and different only in their response to environmental cues. For example, we might expect that many cellular processes are similar across all individuals in the population at least in comparison to, for example, an epithelial versus a B-cell in the tumor microenvironment. Additionally, of course, fungal single cell profiling is in its infancy and better technologies that capture more transcripts per cell will likely produce more distinct clusters. Here, we hope that the reviewer agrees that we have provided evidence that the technology and analyses have succeeded to delineate, at least to some degree, between a handful of different types of cellular responses which we discuss individually below.

Additionally, we would like to clarify the nature of the methodology used here. We use a neural network (scVI) to process the raw data into a lower dimensional latent (hidden space). Leiden clustering is based solely on the expression profiles; the algorithm looks for evidence of separation between groups of cells. UMAP only maps the latent space (from scVI) down to 2D for visualization. Finally, pathway analysis used the identified Leiden clusters; this supervised analysis ultimately yields z-scores. As such, if a signature has higher expression in one Leiden cluster compared to other clusters, pathway analysis used in this manner is a form of cross validation: random clusters are not expected by chance to contain differentially expressed genes from the same functional categories. That is, the clusters point to interesting biological differences.

The microscopy depicted in Figure 2 C-E also represents a second source of cross-validation; it supports the separation of these clusters based only on GFP/RFP strains. It is therefore independent of the sc-transcriptomics data. Of course, this is not definitive and we agree fully with the reviewer that additional statistics can be used to examine stability which we present next.

Author response image 2
Download asset Open asset

We revisited the analysis of UT cells (this corresponds to Results subsection 3, paragraphs 1 and 2 L184-211). The reviewer’s question is primarily regarding the stability of the blue pink and green clusters in the clustering above on the right hand side (Author response image 2, this was Figure 2A in the original manuscript). Towards this end, we labeled every cell with a unique Cell ID (x-axis of the right hand side of Author response image 2). Then, we resampled 100x each time picking a subset of 95% of all UT cells each time, reapplying the scVI algorithm (to re-compute the model of expression) and reapplying the Louvain clustering algorithm (code available at https://github.com/vdumeaux/sc-candida_paper). We then determined which cluster each cell was found in. In the figure on the right hand side (Author response image 2), original cluster 1 corresponds to the darkblue/G1-S phase cluster in the left hand side image, cluster 3 corresponds to the green/M phase cluster and cluster 2 corresponds to sampled cluster labelled `1-3’. Some cells that were placed in cluster 1 originally again find themselves in cluster 1 in the resampling (that is the light purple pattern at the top right). Other cells are placed in cluster `1-3’. This means that, although there is good evidence that cluster 1 exists, the Leiden clustering is far from perfect and there is some confusion with clusters 2. This is also true for original cluster 3. Sometimes these cells are re-classified as cluster 3 but sometimes they are grouped again with cluster `1-3’.So although the three clusters are certainly not absolutely distinct there is strong evidence for all three clusters.

To better address the instability of some clusters within the manuscript, we re-represent our results so that they do not rely on the grouping of cells into discrete clusters. We produced an alternative to Figure 2 (page 8 of new manuscript) which shows the signature scores for each cell individually. For these signatures, we can observe that cells differentially express these signatures and that when a cell expresses a specific signature highly, it also tends to express other signatures highly. For example, a set of cells (cell indices >3000) exhibit elevated expression of genes involved in M phase, the heat-shock response and glycolysis.

Again we assigned all cells a unique Cell ID. They are ordered from left to right in the same manner across all 6 plots, providing the ability to compare them more readily, and to better highlight the differences of expression for specific signatures/biological processes. The text in the subsection describes essentially the three sets of cells (Index 0-2700, 2701-3000, and 3001-5000) and their inter-relationships. Note that here the analysis w.r.t. UT cells no longer depends upon the Leiden clusters.

Perhaps not surprisingly, the Leiden clusters identified post drug treatment (Figure 3, page 11 of new manuscript) are far more stable when we apply the same resampling technique. Figure 3 – Supplement figure 1C has been added to the manuscript. It shows that the 5 main clusters are quite stable (the 19 comet clusters were not included in the results due to their small sizes and given that we considered them both as one pooled comet and as distinct clusters in our analyses).

Figure 3 —figure supplement 1C was included to depict the considerable stability of the clustering that included untreated (UT) and drug-treated cells.

Figure 5B was also added to characterize the stability of clusters identified when investigating FCZ day3 survivors at day 3 and 6 after resuspension in fresh media. Most identified clusters are robust except for cluster 1 and 2 which includes α/Rd cells at day 6 and day 3, respectively. The lack of robust transcriptional differences between cluster 1 and 2 is described L347-350 and confirms that day 6 α/Rd cells likely evolved from day 3 α/Rd cells.

5. I have some concerns about the comparison between untreated and treated at the same (or similar) time points. If the transcript profile changes over time, as presumably it would, then it is reasonable to expect that the treated and untreated cells are at different equivalent times because their growth rates are different. I have a similar concern with comparing different time points (days 3 versus 6) because fresh medium (YPD) was used in this experiment. How much are transcript abundances static over time? This should minimally be discussed further in the manuscript. Also, the conclusions in L472-475 and L579-592 may need this caveat mentioned.

We agree that comparisons between UT and drug-treated populations are problematic. However we do not rely on direct comparisons between UT and drug-treated cells. Rather, analysis of the UT cells was performed to validate the efficacy of the sc-transcriptomics assay and provide a ‘sanity check’.

Unfortunately, the previous version of the manuscript had several sources of confusion with respect to the growth rate of the different C. albicans populations (for example Reviewer #2, Questions 4 and 8), and we have revisited the Discussion points with these criticisms in mind. Briefly, we better emphasize that the primary goal is not to compare UT with drug treated cells and that the major finding of the manuscript is the existence of the ⍺/Rd and β/Sd subpopulations initially identified in FCZ at day 2. This result is independent of any other time point or drug treatment (including untreated). The manuscript focuses on the observation that this same split between ⍺/Rd and β/Sd cells continues to exist at day 3 and even at day 6. While it is possible that some of the distinction between ⍺/Rd and β/Sd cells is not entirely due to drug treatment, it also cannot be due to cell cycle or nutrient levels alone and was not observed with other drug treatments (e.g. CSP and RAPA).

https://doi.org/10.7554/eLife.81406.sa2

Article and author information

Author details

  1. Vanessa Dumeaux

    Department of Anatomy and Cell Biology, Western University, London, Canada
    Contribution
    Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1280-6541

  2. Samira Massahi

    Department of Biology, Concordia University, Montreal, Canada
    Contribution
    Formal analysis, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared

  3. Van Bettauer

    Department of Computer Science and Software Engineering, Concordia University, Montreal, Canada
    Contribution
    Software, Formal analysis, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared

  4. Austin Mottola

    Shmunis School of Biomedical and Cancer Research, The George S. Wise Faculty of Life Sciences, Tel Aviv University, Tel Aviv-Yafo, Israel
    Contribution
    Validation, Investigation, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared

  5. Anna Dukovny

    Shmunis School of Biomedical and Cancer Research, The George S. Wise Faculty of Life Sciences, Tel Aviv University, Tel Aviv-Yafo, Israel
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  6. Sanny Singh Khurdia

    Department of Biology, Concordia University, Montreal, Canada
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  7. Anna Carolina Borges Pereira Costa

    Department of Biology, Concordia University, Montreal, Canada
    Contribution
    Resources
    Competing interests
    No competing interests declared

  8. Raha Parvizi Omran

    Department of Biology, Concordia University, Montreal, Canada
    Contribution
    Resources
    Competing interests
    No competing interests declared

  9. Shawn Simpson

    Department of Computer Science and Software Engineering, Concordia University, Montreal, Canada
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  10. Jinglin Lucy Xie

    Department of Chemical and Systems Biology, Stanford University, Stanford, United States
    Contribution
    Writing – review and editing
    Competing interests
    No competing interests declared

  11. Malcolm Whiteway

    Department of Biology, Concordia University, Montreal, Canada
    Contribution
    Resources, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6619-7983

  12. Judith Berman

    Shmunis School of Biomedical and Cancer Research, The George S. Wise Faculty of Life Sciences, Tel Aviv University, Tel Aviv-Yafo, Israel
    Contribution
    Investigation, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared

  13. Michael T Hallett

    Department of Biochemistry, Western University, London, Canada
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    michael.hallett@uwo.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6738-6786

Funding

Natural Sciences and Engineering Research Council of Canada (RGPIN-2018-05085)

  • Michael T Hallett

Canada Foundation for Innovation (37083)

  • Michael T Hallett

European Research Council (951475)

  • Judith Berman

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank L Pachter and the organizers of BIRS workshop #17w5134 (Oaxaca, Mexico) for several conversations that motivated this effort especially with respect to affordable open biotechnology. We thank CA Jackson and D Gresham for early access to their protocol, C Law for assistance with the microscopy, A Villani for several nice optimizations, and members of the McCarroll lab for their careful advice and outstanding responsiveness.

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Antonis Rokas, Vanderbilt University, United States

Reviewer

  1. David Shore, University of Geneva, Switzerland

Version history

  1. Received: June 25, 2022
  2. Preprint posted: July 21, 2022 (view preprint)
  3. Accepted: October 26, 2023
  4. Accepted Manuscript published: October 27, 2023 (version 1)
  5. Version of Record published: December 6, 2023 (version 2)

Copyright

© 2023, Dumeaux et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 394
    Page views
  • 90
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Vanessa Dumeaux
  2. Samira Massahi
  3. Van Bettauer
  4. Austin Mottola
  5. Anna Dukovny
  6. Sanny Singh Khurdia
  7. Anna Carolina Borges Pereira Costa
  8. Raha Parvizi Omran
  9. Shawn Simpson
  10. Jinglin Lucy Xie
  11. Malcolm Whiteway
  12. Judith Berman
  13. Michael T Hallett
(2023)
Candida albicans exhibits heterogeneous and adaptive cytoprotective responses to antifungal compounds
eLife 12:e81406.
https://doi.org/10.7554/eLife.81406

Categories and tags

Further reading

    1. Computational and Systems Biology

    Homeostasis, injury and recovery dynamics at multiple scales in a self-organizing mouse intestinal crypt

    Louis Gall, Carrie Duckworth ... Carmen Pin
    Research Article

    The maintenance of the functional integrity of the intestinal epithelium requires a tight coordination between cell production, migration and shedding along the crypt-villus axis. Dysregulation of these processes may result in loss of the intestinal barrier and disease. With the aim of generating a more complete and integrated understanding of how the epithelium maintains homeostasis and recovers after injury, we have built a multi-scale agent-based model (ABM) of the mouse intestinal epithelium. We demonstrate that stable, self-organizing behaviour in the crypt emerges from the dynamic interaction of multiple signalling pathways, such as Wnt, Notch, BMP, ZNRF3/RNF43 and YAP-Hippo pathways, which regulate proliferation and differentiation, respond to environmental mechanical cues, form feedback mechanisms and modulate the dynamics of the cell cycle protein network. The model recapitulates the crypt phenotype reported after persistent stem cell ablation and after the inhibition of the CDK1 cycle protein. Moreover, we simulated 5-fluorouracil (5-FU)-induced toxicity at multiple scales starting from DNA and RNA damage, which disrupts the cell cycle, cell signalling, proliferation, differentiation and migration and leads to loss of barrier integrity. During recovery, our in-silico crypt regenerates its structure in a self-organizing, dynamic fashion driven by dedifferentiation and enhanced by negative feedback loops. Thus, the model enables the simulation of xenobiotic-, in particular chemotherapy-, induced mechanisms of intestinal toxicity and epithelial recovery. Overall, we present a systems model able to simulate the disruption of molecular events and its impact across multiple levels of epithelial organization and demonstrate its application to epithelial research and drug development.

    1. Computational and Systems Biology

    Genetics: From mouse to human

    Arya Mani
    Insight

    A deep analysis of multiple genomic datasets reveals which genetic pathways associated with atherosclerosis and coronary artery disease are shared between mice and humans.

    1. Computational and Systems Biology

    Shared and distinct pathways and networks genetically linked to coronary artery disease between human and mouse

    Zeyneb Kurt, Jenny Cheng ... Xia Yang
    Research Article

    Mouse models have been used extensively to study human coronary artery disease (CAD) or atherosclerosis and to test therapeutic targets. However, whether mouse and human share similar genetic factors and pathogenic mechanisms of atherosclerosis has not been thoroughly investigated in a data-driven manner. We conducted a cross-species comparison study to better understand atherosclerosis pathogenesis between species by leveraging multiomics data. Specifically, we compared genetically driven and thus CAD-causal gene networks and pathways, by using human GWAS of CAD from the CARDIoGRAMplusC4D consortium and mouse GWAS of atherosclerosis from the Hybrid Mouse Diversity Panel (HMDP) followed by integration with functional multiomics human (STARNET and GTEx) and mouse (HMDP) databases. We found that mouse and human shared >75% of CAD causal pathways. Based on network topology, we then predicted key regulatory genes for both the shared pathways and species-specific pathways, which were further validated through the use of single cell data and the latest CAD GWAS. In sum, our results should serve as a much-needed guidance for which human CAD-causal pathways can or cannot be further evaluated for novel CAD therapies using mouse models.