Narrow mutational signatures drive acquisition of multidrug resistance in the fungal pathogen Candida glabrata.

Fungal infections are a growing medical concern, in part due to increased resistance to one or multiple antifungal drugs. However, the evolutionary processes underpinning the acquisition of antifungal drug resistance are poorly understood. Here, we used experimental microevolution to study the adaptation of the yeast pathogen Candida glabrata to fluconazole and anidulafungin, two widely used antifungal drugs with different modes of action. Our results show widespread ability of rapid adaptation to one or both drugs. Resistance, including multidrug resistance, is often acquired at moderate fitness costs and mediated by mutations in a limited set of genes that are recurrently and specifically mutated in strains adapted to each of the drugs. Importantly, we uncover a dual role of ERG3 mutations in resistance to anidulafungin and cross-resistance to fluconazole in a subset of anidulafungin-adapted strains. Our results shed light on the mutational paths leading to resistance and cross-resistance to antifungal drugs.

Introduction

Each year, fungal infections affect >1 billion people worldwide and cause 1.5 million deaths. Current challenges to overcome this trend include the lack of fast and accurate diagnoses and the rise of antifungal drug resistance. Acquisition of antifungal resistance is particularly worrying, given the limited number of available compounds. The widespread use of antifungal agents to counteract the high clinical, agricultural, and economic burden caused by various fungal pathogens, coupled with the high ability of fungi to adapt to selective pressures, have resulted in an alarming increase in the rates at which fungal species or isolates resistant to one or multiple drugs are identified., As a result, we are witnessing a global epidemiological change represented by the increased incidence of previously uncommon species with a greater ability to adapt to drugs, the increased failure of therapies due to adaptation of the infecting clone, and the common appearance and rapid spread of deadly outbreaks caused by resistant lineages. These trends affect all major human fungal pathogenic genera, including Candida, Aspergillus, Cryptococcus, and Pneumocystis. Despite the pressing challenge that the emergence of antifungal resistance represents for human health and food security, we have a limited understanding of the evolutionary processes leading to drug adaptation in fungi. Although we know common resistance-conferring mutations and major resistance mechanisms operating in many fungal pathogens, these represent the culmination of an adaptation process. This evolutionary process remains understudied because most of our knowledge derives from already-adapted clones, and from the exploration of a usually limited set of known target genes. In this regard, the use of an in vitro evolution approach coupled to whole-genome sequencing represents a promising research avenue.

Candida species are among the main causes of hospital-acquired fungal infections. C. albicans is the most common cause of candidiasis, but the relative incidence of non-albicans Candida species is on the increase, with C. glabrata often being the second most prevalent cause of infection. C. glabrata belongs to the Nakaseomyces clade and is phylogenetically closer to Saccharomyces cerevisiae than to most other Candida pathogens, which may imply different routes for drug adaptation as compared to other Candida species. Antifungal resistance in C. glabrata is particularly problematic, as this yeast shows a remarkable ability to adapt to both azoles and echinocandins, thus leading to multidrug resistance (MDR).8, 9, 10, 11 Most antifungals commonly used against Candida are azoles (e.g., fluconazole [flz]), fungistatic drugs that inhibit a lanosterol demethylase encoded by ERG11, and echinocandins (e.g., anidulafungin [ani]), which inhibit 1,3-β-d-glucan synthase encoded by FKS genes and are fungicidal to Candida. Most prevalent mechanisms conferring protection against azoles in yeasts involve alterations in the target enzyme or overexpression of drug efflux pumps. Known mechanisms of azole resistance in C. glabrata almost exclusively consist of gain-of-function mutations in PDR1, which encodes a transcriptional regulator of drug efflux pumps, whereas echinocandin resistance has been linked to non-synonymous variations in two conserved hotspot (HS; i.e., frequently mutated) regions of FKS genes. Antifungal drug resistance, tolerance, and adaptation are related to the ability of a cell to respond to stress. Under stress, genome maintenance and repair mechanisms are altered, which may lead to the appearance of resistance phenotypes., Rapid adaptation to varying conditions, including exposure to drugs, has been attributed to a remarkable genomic plasticity in Candida. In C. glabrata, a large degree of genomic and phenotypic variation has been described between and within genetically diverse clades, and even within clonal populations infecting a patient., Previous studies on in vitro-acquired drug resistance in C. glabrata have evaluated the fitness costs of echinocandin resistance or used transcriptomics to unveil the mechanisms contributing to azole resistance, but the genome-wide genetic alterations involved during this process remain elusive. In addition, the genetic underpinnings of MDR in this pathogen are poorly understood.

Results

C. glabrata has a widespread ability to acquire drug and MDR

Here, we set out to explore the evolutionary adaptation of C. glabrata to azoles and echinocandins using an in vitro evolution approach coupled to phenotyping and targeted gene and whole-genome sequencing (Figure 1; STAR Methods). To this end, 12 strains representing the 7 previously described C. glabrata clades were subjected to increasing concentrations of antifungal drug(s) in the following regimes: fluconazole (FLZ samples; note the use of uppercase letters for samples/conditions as opposed to lowercase letters for the drug), anidulafungin (ANI), and both drugs in combination (ANIFLZ). In addition, to gain insight into mechanisms of cross-resistance, adaptation to serial exposure to both drugs was studied by growing isolates from the final steps of the ANI samples under the flz regime (AinF) and, conversely, final FLZ isolates under ani (FinA). Finally, control populations of all of the strains were grown for the same time without any drug (YPD). The experiment comprised a total of 288 independently evolved populations. When exposed to a single drug or to the two drugs in a sequential manner, all of the populations survived the entire experiment. However, when simultaneously exposed to both drugs, 21 populations (43.75%) died, including all replicates of each of 2 strains from clade I (CST109) and clade III (M12). Nevertheless, populations from other strains from these clades survived, indicating that low adaptation potential is strain- and not clade-specific. We analyzed available parental sequences of the two strains unable to adapt to ANIFLZ and found that they shared eight genes (the S. cerevisiae orthologs of SWI6, CDC3, LAP2, MAD1, MNN4, RSN1, and SQS1 and the gene CAGL0C05313 g) with alterations that were not present in the parentals of the surviving strains within the same clades (Table S1).

Figure 1

Schematic representation of the in vitro evolution experiment

A total of 48 populations, quadruplicates of each of the 12 strains, were grown with increasing concentrations of flz (FLZ samples), ani (ANI), both drugs in combination (ANIFLZ), and no drug (YPD). Subsequently, ANI samples were grown in flz (AinF), whereas FLZ samples were grown in ani (FinA). The experiment involved batch serial transfer of the samples every 3 days, in which every second passage involved an increase in drug concentrations up to 4 and 196 μg/mL ani and flz, respectively (Table S4; STAR Methods). After the final passage, an aliquot was plated for single colony isolation and storage.

We determined susceptibility using the minimum inhibitory concentration (MIC) and the relative area under the curve (rAUC) measurements (Figures 2A, 2B, S1A, and S1B; STAR Methods). All of the surviving strains acquired stable resistance to the exposed agent(s); that is, the resistance phenotype was kept for several generations in standard growth conditions after the removal of the selective agent (Figures 2B and S1C; Data S1), indicating that the phenotype is genetically encoded. Unexpectedly, we observed increased resistance to flz in a large subset of ANI samples (21/47, MIC > 256), thus showing that adaptation to ani can frequently induce cross-resistance to flz. The reverse process, acquisition of resistance to ani in FLZ samples, was not observed (Figure S1C). Increased resistance to both drugs (MDR) was often achieved, including all surviving ANIFLZ samples, a majority of AinF (91.6%) and FinA (97.9%) samples, and, due to the mentioned cross-resistance, in 44.7% of ANI samples (Figure 2C; Data S1). In serial drug-exposure experiments, previously acquired resistance was rarely lost during exposure to the second drug (1 FinA and 4 AinF samples), indicating that the phenotype is stable. To assess cross-resistance to other antifungal drugs, we tested the growth of a selected panel of evolved strains on other antifungal drugs (Figure S2D). Similar results were observed for the two tested echinocandins (ani and caspofungin), while the two tested azoles presented more disparate patterns, with few strains growing better on voriconazole (vrz) than on flz (discussed below). None of the tested strains presented improved growth on flucytosine (5-FC, pyrimidine analog) or amphotericin B (ampB, polyene) when compared to wild-type (WT) strains, although a few strains presented higher susceptibility to ampB. We evaluated the fitness costs of acquired resistance using AUC values of growth curves in the absence of the drug as a proxy for fitness (fAUC) relative to the fitness of the unevolved (WT) strain (Figure 2D; Data S1). All flz-exposed samples showed a tendency to reduce fitness (p < 10−5, Kolmogorov-Smirnov test), while the mean fitness of ANI samples remained unaltered (p > 0.05). Consistently, a small but significant negative correlation between resistance (rAUC) and fitness levels for flz, but not for ani, was detected (Figure 2E). Nevertheless, many of the flz-exposed samples retained fitness levels within 2 standard deviations of the mean of YPD-exposed strains (56% of ANIFLZ, 77% AinF, 81% FLZ, and 68% FinA), and only a few samples (2.9%, 5/8 of them ANIFLZ) had severely reduced fitness levels below 50% of the corresponding WT strain. These results indicate that resistance, including MDR, is often achieved at mild fitness costs. Finally, we evaluated the repeatability of the fitness and susceptibility outcomes in the parallel evolution experiments for replicates and strains subjected to similar conditions. We did so by comparing the distribution of pairwise differences between samples with respect to assayed fitness and susceptibility levels. Our analysis (Figure S1E) indicates that repeatability may be unique to each phenotype and condition, where AinF and ANIFLZ samples have particularly higher phenotypic variability. In addition, we found that variability was similar among evolved samples of the same or different strains (Figure S1E), suggesting that different strains reached similar phenotypes. Interestingly, we found some exceptions (including the fitness and flz resistance in YPD-evolved and the fitness of FinA-evolved samples) in which the evolved phenotypes are more consistent among samples of the same strain (Figure S1E).

Figure 2

Fitness and drug resistance

(A) We measured relative fitness (the ratio between fitness in each drug concentration versus the no-drug condition [control]) in a time course experiment at several concentrations of flz and ani. Fitness was measured as the area under the time-versus-optical density (OD) curve (fAUC). The graph depicts an illustrative example of two independently evolved replicates of the CST109 strain in the ANI and YPD evolution experiments. The shaded areas represent the median absolute deviation across technical replicates. As a proxy for drug resistance, we defined rAUC as the AUC of these data (normalized by the maximum AUC, in which fitness is maintained across all the range of concentrations [AUCMAX]).50% of growth inhibition, as compared to the no-drug control, is marked as MIC50.

(B) rAUC for flz (top) and ani (bottom) across all samples in our experiments. Each point corresponds to an independently evolved biological replicate. Note that some samples have an rAUC above 1.0, where fitness did not drop upon increasing drug concentration (suggesting high resistance). In addition, Figure S6 includes information about the drug resistance levels among samples with different mutations.

(C) The relationship between ani and flz resistance across all samples. Dashed lines indicate median rAUCs levels for each drug in the YPD samples and rAUCMAX (1.0). Each point corresponds to a biological replicate, and the error bars reflect the median absolute deviation across technical replicates. Each marker corresponds to a different strain.

(D) Fitness in the absence of drug (measured as the log2 fold change in fAUC (see [A] between each sample and the median fAUC in the WT of the matching strain). Note that Figure S6 includes information about relative fitness levels among samples with different mutations.

(E) Fitness in the absence of drugs is slightly correlated with the levels of flz, but not ani, resistance (rAUC). Spearman’s correlation coefficient (r) and p value are shown for flz (left) and ani (right) resistance. The correlation for flz resistance was maintained when considering only samples with mild fitness defects (fitness >−1, r = −0.22, p = 0.0029). Only resistant samples, defined as those with a log2 fold increase above 1 as compared to the WT (Figure S1D), were included in this analysis. The individual fitness and susceptibility measurements for each sample can be found in Data S1.

The FKS mutational spectrum in resistant strains expands beyond HS regions

We used a target sequencing approach to screen 121 ani-adapted strains for mutations in the typically surveyed HS of FKS genes (Data S2; STAR Methods). In addition, we selected 77 representative (considering clades, susceptibility levels, and FKS mutations) samples for whole-genome sequencing and called small variants (SVs), copy-number variations (CNVs), and genomic rearrangements (GRs) appearing de novo in each of the evolved samples (Data S3; STAR Methods). All 121 ani-evolved strains presented newly acquired non-synonymous (ns) mutations in the targeted FKS regions (Data S2), which indicates that FKS mutations may be necessary for ani adaptation. Mutations preferentially occurred in FKS2 over FKS1 and in HS1 over HS2 (Figure 3), suggesting a more prevalent role of these loci. Notably, 22% of FKS mutations were outside the HS regions. Three resistant strains carried only such non-HS FKS mutations (FKS1-R1422L and FKS1-F708S; FKS1-W681L and FKS2-K265∗; and FKS2-A651T; Data S2), and whole-genome sequencing of these strains revealed no additional mutations outside FKS genes that could explain their resistance phenotype (see below). These observations suggest that some of these non-HS FKS mutations contribute to resistance and emphasize the importance of studying FKS genes beyond HS regions. In addition, we tested whether the distance of non-HS mutations to the actual HS is related to the level of ani resistance in samples harboring only non-HS mutations. We could not find any such significant correlation (Spearman rho = −0.14, p < 0.11), suggesting that non-HS and HS mutations confer similar levels of resistance. Overall, the most frequently mutated site in ani-adapted samples was FKS2-F659 (63 samples, 52.1%; Data S2), with the most prevalent alteration being F659del (52 samples, 43%), which was the only FKS mutation in 26 samples (21.5%). This finding suggests that, as compared to replacements, amino acid deletions may more efficiently prevent the binding of the drug, and reinforces the need to consider this type of mutation. Finally, 26 samples exposed to ani (19.8%) carried a truncation in one of the FKS genes (2 of them with a GR breaking the coding region (Figure S3; STAR Methods) in combination with a ns mutation in the other paralog, indicating that this specific combination may facilitate adaptation.

Figure 3

Mutational analysis of FKS regions

(A) Distribution of the mutations in studied regions of FKS. A non-negligible presence of mutations outside HSs can be observed. Note that Table S5 includes the oligos used for the sequencing. In addition, Data S2 includes the precise mutations.

(B) Distributions of samples according to the presence of mutations in particular FKS gene and distribution of samples according to the presence of mutations in FKS HSs.

(C) Mutational signatures per sequenced regions: FKS1 and FKS2_1 and FKS2_2. Mutated positions are shown as highlighted boxes at the corresponding amino acid in the mutation, over a gray background. Color scale, from white to red, indicates the observed number of mutations (log scale). Darker gray boxes indicate HSs and the white-framed box in FKS2_1 marked positions for other possible mutational HSs. The bottom part of the graph represents an enlargement in HSs and mutations in their close proximity.

Mutational landscapes in resistant strains reveal a high diversity of genetic alterations affecting a restricted set of recurrently mutated genes

The analysis of genome-wide mutational patterns revealed no newly acquired SVs in YPD samples, while the drug-evolved strains accumulated a small number (<10) of variants (Figures 4A and 4B). This indicates that susceptible strains are a few mutational steps away from acquiring resistance. Strains carrying distinct MSH2 variants (Figures 4A and 4B) did not accumulate a different number of mutations, thereby supporting the notion that these represent natural, functional variants rather than hypermutator mutations. As expected,, we found that aneuploidies were common in experiments involving exposure to flz, but they were not detected in cells exposed only to ani (Figure 5A). Total or partial aneuploidies in chromosome E (ChrE), encompassing ERG11, were the most common, appearing in 11/16 FLZ, 4/15 AinF, and 2/6 ANIFLZ samples. Most (10/11) FLZ samples with the ChrE aneuploidy retained it upon further exposure to ani (FinA). One strain presented a partial ChrE aneuploidy resulting from unbalanced translocation with ChrJ (Figure S3D; STAR Methods), suggesting that GRs can drive drug resistance. Importantly, we detected no heterozygous variant in any of the duplicated chromosomes, indicating they have not accumulated new mutations since their duplication, and, therefore, that aneuploidies were adaptive per se and not because they allowed faster evolution of duplicated genes. To investigate whether aneuploidies conferring flz resistance were rendering strains avirulent, we used an in vivo Galleria mellonella model (STAR Methods) to assess the virulence of a WT strain and two of its descendant FLZ strains, one of which presented chromosomal duplications in ChrE and ChrI. Our results (Figure 5B) show that all of the descendant strains remained virulent, suggesting that flz resistance or the presence of aneuploidies are compatible with virulence.

Figure 4

The number of small variants (synonymous and non-synonymous) that appear during the experiment

(A) To select only newly acquired mutations in each drug-evolved sample, we subtracted from called variants those also called in the corresponding WT, YPD, and the parental drug condition (ANI for AinF, and FLZ for FinA), while the corresponding variants called in WT, ANI, AinF, FinA, and FLZ samples were subtracted from those found in the YPD sample. The dashed lines, from bottom to top, correspond to 1 and 5 mutations, respectively. We also represent the presence of ≥1 ns variants in the MSH2 gene in the WT strain. The bars represent the mean number of mutations across biological replicates and the error bars represent the standard deviation.

(B) The same as in (A), but showing the fraction of protein-altering mutations.

Figure 5

The role of aneuploidies in drug resistance

(A) We calculated the median relative coverage per gene for all samples analyzed in this work. This parameter appeared to be correlated with the distance to the telomere (STAR Methods), so that the log2 ratio to the YPD (of the corresponding strain) was used as a proxy for the gene copy number. Shown is the rolling-median of this value for windows of 50 genes and chromosomes where large duplications were observed (chromosomes E, I, A, and L). Data for chromosomes I, A and L are shown only for those strains in which aneuploidies are observed. Each column corresponds to a sample (ordered as in Figure 6), and the “∗” and “X” correspond to FinA samples in which the parent (FLZ) aneuploidy was maintained or lost, respectively. ERG11, PDR1, and TPO3 are genes that we speculate could be driving the selective advantage of the aneuploidy (see Results). All of the values were cut off at 1.0 (2× coverage as compared to the YPD) for clarity.

(B) Survival of Galleria mellonella larvae during 6 days after inoculation of EB0911 (WT strain) and 2 flz resistant progenies: 3B_FLZ (without aneuploidies) and 3H_FLZ (presenting both ChrE and ChrI duplications).

Figure 6

Aneuploidies and recurrently mutated genes

Each drug is associated with a particular set of mutated genes and aneuploidies. Columns represent the evolved samples, each strain indicated by a number: 2, CST34; 3, EB091; 4, CST78; 5, M12; 6, EF1237; 7, EF1620; 8, F15; 9, CBS138; 10, P35; 11, BG2. Replicates of the same strain appear in the same order as in the experimental plate. Colors indicate the experimental condition. Blocks show, from top to bottom, chromosomal alterations, mutated genes, and susceptibility data. Whole and partial (P) chromosomal duplications appearing newly in each condition are marked as red, while losses are marked as light salmon boxes. Protein-altering mutations (gray boxes) and losses (black boxes) of genes appearing in at least 2 drug-evolved samples are shown. Note that we found a balanced translocation in FKS1 (T) and a deletion in the ERG3 promoter region (Pr) (Figure S3; Results; STAR Methods). PTC stands for premature termination codon. Pink arrows indicate the parent-daughter relationships for 3 AinF samples that did not present any new alteration in recurrent genes. Note that Figures S3 and S4 and Data S3 provide additional information about these mutations and genomic rearrangements. In addition, Figure S6 shows the association between these mutations and fitness or drug-resistance levels.

To identify mutations likely associated with the resistance trait, we selected genes that were mutated at least twice independently in our experiment. This search identified nine genes (ERG11, PDR1, CDR1, CNE1, EPA13, FKS1, FKS2, ERG3, ERG4; Figure 6). Importantly, all of the resistant strains carried mutations or duplications in at least one of these genes, and the subset of mutated genes largely separated samples by treatment. This strong association of acquired mutations, treatment, and phenotypes indicates that a limited set of genes is central for the acquisition of resistance. The most common altered gene under exposure to flz was PDR1, which was in many instances (14/37 strains) accompanied by alterations in ERG11 (Figure 6; Data S3). Although less common, five resistant strains contained no PDR1-related mutations or aneuploidies (Figure 6), indicating that alternative mechanisms confer resistance on their own. These strains harbored mutations in ERG3 (3 AinF strains, discussed below) and ERG11 (2 strains). Importantly, ERG11 mutations and aneuploidies in ChrE, bearing this gene, were strongly anticorrelated, with a single ANIFLZ sample carrying both alterations. In this case, the mutation was present in the two alleles, suggesting that the mutation preceded the chromosomal duplication. Among ERG11 mutations, K152 was the most altered amino acid (12/16 samples), followed by ERG11-Y141 (2/16 samples). Although common in other Candida species, these mutations have not been commonly reported in C. glabrata. Structural analysis revealed that both altered residues were close to the azole binding pocket (Figure S4).

We next assessed whether the catalog of mutations found in our in vitro analysis was representative of what can be found in clinical strains. To this end, we compared this catalog with variants found in 393 C. glabrata clinical isolates with genomes publicly available at Candidamine (https://candidamine.org/). Our results (Figure S5; STAR Methods) show that the overlap of specific mutations is very low. This low overlap is, however, expected from the actual large diversity of the identified mutations in our experiments (Figure S5B; Data S2 and S3) and is similarly low for mutations identified in actual clinical surveys (e.g., SENTRY). These results suggest that although the set of genes recurrently mutated during the acquisition of resistance is rather limited (nine genes), the number of specific mutations (i.e., which residue is mutated and what type of mutation occurs) is large and highly diverse and only partially covered by our experiment or clinical surveys.

Decreased fitness of some resistance-conferring mutations could hamper their detection in the clinics, as clinical isolates are not obtained in selective conditions. To explore possible fitness trade-offs of specific mutations, we evaluated whether strains harboring each type of mutation had a particular fitness or susceptibility level. Consistent with the fitness results presented above, most of the mutations had no significant effect on fitness in the absence of the drug (Figure S6). However, we found that strains harboring CNE1 truncations or ChrL and ChrA duplications presented lower fitness, indicating that some resistance mechanisms may generate decreased growth (Figure S6). On another note, we found that most strains had similar flz and ani susceptibility levels independently of the mutation type (i.e., we found no differences in flz resistance among strains with ERG11 mutations or ChrE duplications) (Figure S6). Finally, we investigated whether there was a correlation between the number of different genes with acquired mutations and fitness/susceptibility levels in any of the evolution conditions. We found no significant Spearman correlation (p < 0.05) after removing a single outlier AinF sample with a particularly high number of new mutations and low ani resistance. These results indicate the lack of a general correlation between the numbers of acquired mutations and these evolved phenotypes. Our data suggest that different evolutionary paths drive similar levels of drug resistance and fitness in a strain-independent and mutation-independent manner.

Crosstalk between echinocandin and fluconazole resistance

In the experiments of sequential exposure to the two drugs, all of the samples successfully adapted, in turn, to the two challenges. When adapting to the new drug, most samples (90 of 95) retained the previously acquired resistance, resulting in MDR (Figures 7A, 7B, S1C, and S1D). However, three sequenced samples lost the previously acquired resistance upon the change in selective conditions (according to MIC, see Figures S1C and S1D). These included a FinA sample and two AinF samples. This FinA sample acquired a premature termination codon in PDR1, which may revert the flz resistance conferred by previous mutations in this gene. In the two AinF samples that lost resistance to ani, we found frameshift mutations in FKS2 downstream of the ani resistance-conferring mutations inherited from the parental ANI samples (Figure 7A). Interestingly, both of the ANI parents carried only one FKS2 mutation and alterations in CNE1 ortholog, encoding an endoplasmic reticulum (ER) protein involved in the quality control of misfolded proteins. This remarkable coincidence suggests that the combination of these alterations results in a higher propensity to lose resistance, although this hypothesis needs further study. Except for a single ChrA duplication found in one strain, most ANIFLZ samples showed mutational signatures similar to those acquired during sequential exposure to the two drugs (AinF and FinA; Figure 6). This observation suggests that the genetic basis driving the acquisition of resistance to each of the drugs is similar when the two drugs are in combination.

Figure 7

ERG3 mutations and multidrug resistance

(A) Biplot showing the relationship between resistance (rAUC) toward ani and flz for a series of ANI/AinF related samples. The gray dashed lines indicate the rAUC = 1.0 (where fitness is maintained across the range of concentrations; Figure 1A) and the median rAUC across YPD samples for each of the 2 drugs. Each sample is represented by a symbol, with the color indicating the sample type: ANI (pink) and AinF (red) samples. The pink dashed lines indicate parent-daughter relationships (ANI-AinF) between the samples. The symbols represent different types of ERG3 mutations, and the gray circles outline 3 samples that did not acquire any new mutation in the recurrent genes in AinF. The 2 ANI samples with alterations in CNE1, which lost ani resistance due to truncations in FKS2(∗) in AinF samples, are marked. One of the ANI samples showed high ani resistance (above 1.0, meaning the fitness was higher in ani than in no drug), but also showed low basal fitness, which means that the high resistance value may be not representative. Error bars reflect the median absolute deviation across technical replicates.

(B) Relationship between rAUC of ani and flz in FLZ (light blue) and FinA (dark blue) samples. The green dashed lines indicate parent-daughter relationship (FLZ-FinA). The gray dashed lines indicate the rAUC = 1.0 (where fitness is maintained across all the range of concentrations; Figure 1A) and the median rAUC across YPD samples for each of the 2 drugs. No acquisition of ani resistance was observed in FLZ samples but only as a result of ani (FinA). The symbols represent the presence of ERG11 missense mutations or chromosome E aneuploidies. Two FinA samples showed a drop in flz resistance levels. One of them carried a PDR1 premature termination codon (∗), which resulted in susceptibility according to our MIC-based thresholding (STAR Methods) and reduced flz resistance below the median rAUC value of YPD samples. The other sample carried ERG4 mutation that resulted in a reduction but not a total loss of flz resistance. Error bars reflect the median absolute deviation across technical replicates.

(C) Non-synonymous (including missense and STOP loss) ERG3 mutations are associated with higher flz resistance (rAUC) in ANI samples. The p value corresponds to a Kolmogorov-Smirnov test. The corresponding AinF and FLZ samples are also shown for comparison of flz-resistance levels. The dashed symbols represent samples that were found to be flz susceptible according to our MIC-based thresholding (STAR Methods). Note that 2 samples (marked with “?”) were found as susceptible but have rAUC values in the range of resistant samples. This mismatch is clarified in Figures S7C and S7D. In addition, see Tables S2 and S3 for further information on the ERG3 mutations found in each sample.

(D) The presence of ERG3 non-synonymous mutations is correlated with discrete flz resistance in ANI samples. The number of ANI samples in each category and the p value of a Fisher test are shown.

(E and F) Growth competition between ani-resistant strains with and without ERG3 mutation (note that Table S5 includes the oligos used for sequencing). The y axis presents the calculated ratio of a sample with mutated ERG3 gene and the x axis ratios aimed at the beginning of the experiments. The error bars represent the standard deviation across technical replicates. (E) In vitro fitness competition of 2 pairs of strains: 1-CRISPR transformant ERG3 (D122Y) versus CRISPR transformant ERG3(WT) with NAT1 and 2-CRISPR transformant ERG3 (D122Y) with NAT1 versus 3H_ANI (ERG3 WT). The competition was conducted over a 24-h period and in YPD and YPD supplemented with 0.5 μg/mL ani in triplicates. (F) Two independent competition experiments in vivo. The fungal burden was obtained from 3 separate larvae for each of the initial mix of populations.

A remarkable finding of our experiment is the cross-resistance to flz found in a significant fraction of ANI samples (see above). Whole-genome sequencing of 7 of these strains revealed that all of them carried alterations in ERG3, which encodes the C-5 sterol desaturase of the ergosterol biosynthesis pathway (Figure 6). This association was further explored by Sanger-based target sequencing of the ERG3 gene in the remaining ani-evolved strains, which showed that all 21 ani-evolved strains showing cross-resistance to flz (MIC > 256 μg/mL) carried alterations in ERG3 (Table S2). Accordingly, we detected a significant association between ERG3 ns mutations and flz resistance in ANI samples (Figures 7C, 7D, and S7A–S7D). Interestingly, these samples showed lower levels of flz resistance when compared to FLZ samples (Figure S7C). This finding indicates that the quantitative contribution of ERG3 mutations to flz resistance differs from that of PDR1 or ERG11 alterations and suggests different mechanisms of resistance in FLZ and ANI samples. When ERG3-mutated strains were subsequently exposed to flz (AinF), three of them did not acquire additional mutations in PDR1 or ERG11, nor did they present ChrE duplications, thereby suggesting that their ERG3 mutations were sufficient for their survival in flz. In support of this notion, the levels of flz resistance of these three AinF samples and their respective ANI parents were similar (Figure 7A). However, the relationship between ERG3 alterations and cross-resistance to flz was incomplete and mutation dependent. We found that of 28 ANI samples harboring ERG3 mutations, 6 carrying premature stop (3), missense (2), and frameshift (1) mutations retained WT levels of susceptibility. The absence of resistance in strains carrying ERG3 mutations leading to truncated proteins is compatible with earlier work showing that ERG3 deletion in C. glabrata does not affect flz susceptibility. Consistent with some ERG3 alterations being selected under exposure to ani, 2 ANIFLZ and 6 FinA samples bearing ERG3 changes additional to PDR1 and/or ERG11 mutations were detected (Figure 6). Incidentally, another FinA sample carried a deletion in the gene immediately upstream of ERG3 (CAGL0F01815 g, of unknown function), which we speculate may result in regulatory alterations of ERG3 through disruption of the promoter (Figure S3A; STAR Methods). To investigate the relationship between ERG3 mutations and flz resistance further, we re-introduced one of the ERG3 mutations (D122Y) into two WT strains (CBS138 and EB0911) and an ani-evolved and flz-susceptible progeny of EB0911-3H_ANI. In addition, we reverted ERG3 to the WT sequence in one strain (3B_ANI, progeny of EB0911) originally harboring ERG3 D122Y mutation. We then assayed the susceptibility phenotype of these transformants and the original strains. Our results (Figures S2A and S2B) show that the introduction of the D122Y mutation in ERG3 led to increased resistance to flz, and that the reversion of the mutation had the opposite effect, confirming the link of ERG3 and flz susceptibility. We noted that the effect of this mutation was stronger in an ani-resistant background as compared to a WT background, where growth on flz was observed at a later time point. Our results support a dual role of ERG3 alterations in the adaptation to ani and in the development of cross-resistance to flz in C. glabrata.

To gain mechanistic insight into these ERG3 alterations, we performed various experiments. We tested whether the introduced ERG3 alterations were associated with altered response to various stresses. Our results (Figure S2C) suggest no major effects, with the exception of a lower tolerance to membrane (SDS) and oxidative (H2O2) stresses restricted to a particular ani-resistant background strain. In addition, we traced the order of appearance of ERG3 and FKS mutations along intermediate generations in ANI strains and found equal numbers of cases (2 each) in which either ERG3 or FKS mutations predated the other one, and 5 cases in which both mutations are traced to the same intermediate generation (Table S3). These data suggest that one mutation does not necessarily predate the other one. Resistance to flz is often spontaneously acquired in C. glabrata by partial or total loss of mitochondrial DNA, rendering a so-called petite phenotype. However, we can discard this effect in the identified ERG3 mutants due to the absence of deletions in the mtDNA (Figure 6; Data S3) and the absence of a petite phenotype (Figure S7E). We further analyzed competitive fitness between ani-resistant strains with and without ERG3 mutations using in vitro and in vivo (G. mellonella) competition assays (STAR Methods). Both assays provided similar results (Figures 7E and 7F), supporting a competitive disadvantage of the ERG3 mutants in the absence of the drug. However, when the in vitro competition experiment was performed in the presence of ani, the ERG3 mutant outcompeted the WT. These results support the selection for ERG3 mutations only during ani treatment and point to a possible explanation for the lack of clinical cases showing this alteration.

Discussion

Our study adds support to the suitability of in vitro approaches to study the evolutionary acquisition of resistance to antifungal drugs,,, and contributes to a better understanding of the mechanisms of drug adaptation in C. glabrata. Given the high number of replicates and the drug-specific patterns we consistently observed, we can conclude that the discussed mutations are likely related to the specific drug exposure and not to the experimental setting. Our results show that C. glabrata exhibits a remarkable capacity to acquire resistance to the tested drugs, independently of the phylogenetic background of the strain. This is also true for the case of serial exposure to the two drugs, to which all strains and replicates adapted. However, the combined exposure to both drugs prevented adaptation in a significant fraction of the cases, with two strains from two different clades showing an inability to develop resistance in this scenario. Our results show that neither phylogenetic clade nor the presence of non-synonymous mutations in MSH2 are good predictors of the ability to develop MDR, which is pervasive in C. glabrata. Whole-genome sequencing revealed a relatively limited catalog of a few genes that are commonly affected upon sustained adaptation to antifungal drugs. We observed the appearance of commonly reported alterations in FKS, PDR1, and ERG11 genes, which indicates that our experiment reflects processes that also occur in the clinics. However, 5 other genes (CDR1, CNE1, EPA13, ERG3, and ERG4) were recurrently mutated in our experiments. This finding indicates that alternative mechanisms may be concomitantly used to achieve a stable resistant phenotype. Alterations in the promoter region of the efflux pump CDR1 have already been reported in azole-resistant strains,, and our results suggest that alterations of the protein product may also contribute to flz adaptation. We propose that the observed CDR1 mutations increase azole efflux and thus decrease flz susceptibility. As discussed, CNE1 is involved in the quality control of misfolded proteins in the ER. EPA13 is a sub-telomerically encoded lectin-like adhesin with a role in cell adhesion, whose potential role in drug resistance is unknown. Altered adhesion has been linked to azole resistance in C. glabrata, which may explain why EPA13 deletions could be adaptive under exposure to both azoles and echinocandins in our experiments. ERG4 is another gene involved in the ergosterol biosynthesis pathway, which, similar to ERG3, may influence resistance to flz. Future experiments should help determine the order of appearance of these mutations and their specific roles in drug resistance or adaptation. In addition, our results suggest that GRs and CNVs around these genes are related to drug resistance, as previously proposed in C. albicans. This indicates that the traditional focus on SNPs is underpowered to understand the genomic drivers of drug resistance. Finally, our results suggest that although the set of genes altered during the process of adaptation may be limited, the diversity of possible resistance-conferring mutations in each of the affected genes is very large.

An important result from our experiment is the observation that adaptation to ani often results in cross-resistance to flz (but not the other way around). This result was unexpected, given the different modes of action of the two drugs, in which ani affects the cell wall in a fungicidal manner and flz affects the cell membrane, causing growth arrest. This observation is of high relevance, given the expanding MDR in C. glabrata and also considering that some recent guidelines (e.g., from the Infectious Disease Society of America) recommend an echinocandin-based initial therapy against most invasive Candida spp. infections. Importantly, these findings are consistent with a recent report of flz cross-resistance in ani-adapted C. glabrata isolates. We consider that our results can inform future clinical trials or therapy guidelines. For instance, our data suggest that flz resistance may be common after the failure of ani therapy, so that flz treatment following ani may also result in therapy failure. Thus, monitoring of flz resistance after ani therapy, or the use of a different drug as a second line of therapy may be recommended. Similarly, our results point to the absence of cross-resistance to ani when flz is used as a first therapy or to a high clearance potential of the concomitant use of flz and ani, which may be considered in specific cases. Importantly, many flz-resistant strains were susceptible to vrz, which could be a promising therapeutic alternative. However, this observation was drawn from a few samples and requires further research. The scarcity of sequenced genomes for MDR clinical strains and the lack of information of the treatment regime they were exposed to (STAR Methods) prevented us from assessing how commonly this cross-resistance mechanism occurs in the clinics, something that deserves further investigation. We studied the possible molecular basis of such cross-resistance and found compelling evidence of the involvement of ERG3 mutations. In our experiment, alterations in this gene often appeared under ani exposure and were retained in subsequent flz exposure, sometimes without any further mutation being acquired that would explain the acquisition of resistance to flz. In addition, ERG3 mutations were always present in ani-evolved strains that showed cross-resistance to flz, and we confirmed the causative association of flz resistance of the ERG3 alteration by reintroducing it in a flz-sensitive background. Competition assays between strains carrying the WT and the mutated ERG3 allele showed a competitive disadvantage of ERG3 mutants in the absence of drug treatment, but an advantage in the presence of ani. This underscores the complex fitness trade-offs of resistance-conferring mutations and suggests that the frequency of resistance-conferring alleles is likely to fluctuate after treatment. An intriguing possibility is that clones carrying resistance-conferring mutations and causing therapy failure may be missed during the process of strain identification, as blood cultures and colony isolation is generally performed in the absence of drug exposure. Such phenomenon could partly explain the observed discrepancies between resistance levels of clinical isolates and therapy failure.

Importantly, the link between ERG3 and cross-resistance may not be restricted to C. glabrata as ERG3 mutations leading to the depletion of ergosterol and the accumulation of less toxic sterols when ERG11 is inhibited have been implicated in cross-resistance between azoles and polyenes in S. cerevisiae and C. albicans40, 41, 42, 43 and between echinocandins and azoles in C. parapsilosis., In addition, acquisition of ERG3 mutations upon echinocandin exposure has also been described in C. auris. Why ERG3 mutations are often acquired under exposure to ani and how they contribute to resistance to flz remain unclear and need further attention. A speculative scenario is that certain ERG3 mutations lead to alterations in the membrane composition in a way that partially compensates cell-wall alterations induced by ani exposure. In this regard, it has been reported that cell membrane modifications related to changes in ergosterol production affect the structure and composition of the cell wall.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Toni Gabaldon (toni.gabaldon@bsc.es).

Materials availability

Material generated in this study is available upon request from the lead contact.

Experimental model and subject details

Microbe strains

The 12 strains of C. glabrata used in this study are listed in the key resources table. Eleven clinical strains had been previously analyzed for several phenotypic characteristics, including susceptibility to various drugs. In addition, they have been shown to belong to seven genetically distinct clades. The remaining strain (SLL2_glab) was isolated from an oral wash of a healthy individual from Spain, and can thus be considered commensal. SLL2 glab was sequenced within this project and assigned to clade 7.

Galleria mellonella

Unsexed Galleria mellonella larvae were purchased from DNAT ecosistemas (https://www.dnatecosistemas.es).

Method details

In vitro evolution

We conducted experimental evolution experiments using a batch serial transfer approach (Figure 1). Wild-type (WT) strains were collected from glycerol stocks, plated, left to grow until single colonies could be detected and re-plated again for an overnight culture (YPD agar plate at 37°C). A few colonies were suspended in sterile water and diluted to 2.5 × 105 colony forming units per milliliter (CFU/mL). A 96 deep-well plate (2.2mL) with 450 μL of YPD – the master plate - was inoculated with 50 μL of the cell suspension in four replicates for each strain. To ensure lack of cross contamination the inoculations were organized using a checkerboard design (Figure 1) and visually inspected for unwanted growth in non-inoculated wells. Each well of the deep well plate also included a glass bead to ensure proper oxygen transfer and prevent the samples from sedimentation. The master plate was covered with a sandwich cover (Enzyscreeen BV) to ensure optimal oxygenation and limit evaporation. It was then shaken at 300 rpm, and incubated at 37°C for 72 h. Afterward, 50 μL of each culture was transferred to a fresh 450 μL of YPD medium and left again to grow in the same conditions. Next, 50 μL of samples from the master plate were distributed into four independent 96-well plates containing 450 μL of YPD medium supplemented with the following: 1) an echinocandin: anidulafungin (drug: ani, outcome samples: ANI); 2) an azole: fluconazole (flz, FLZ); 3) anidulafungin and fluconazole (aniflz, ANIFLZ); or 4) no drug (YPD). Adaptation to the drugs involved passages of the (50 μL) samples to a fresh (450 μL) medium every 3 days, and in every second passage the concentrations of flz and ani were gradually increased from 4 μg/mL and 0.016 μg/mL to 192 μg/mL and 4 μg/mL, respectively (Table S4), except YPD where no change in the composition of the medium was applied. For each passage the medium with antifungals was freshly made on the same day using a frozen stock of the drugs. Before each increase in drug concentration, part of the culture was frozen and stored at −80°C (100 μL of the sample in 100 μL of 50% glycerol). All in all, the experiment involved 6 days of adaptation to the same conditions before increasing the stress, and further adaptation. Starting with 4 μg/mL flz and 0.016 μg/mL ani, the experiments finished after 54 days, 18 passages with drugs, and 9 increments in drug concentrations. We estimate this period to involve between 60 to 500 generations (assuming a minimum of three doublings per passage in a 1:10 dilution and a maximum of 5-10 generations/day based on earlier studies). From the last passage we selected, stored and analyzed single colonies that were picked from agar plates and regrown on liquid medium supplemented with the last concentrations of the drugs used in each condition. In the second part of the experiment, we repeated the evolution experiment, this time evolving ANI isolates in flz (AinF), and FLZ isolates in ani (FinA), using the same regimes as explained above. Due to the inability to re-grow two samples (1 ANI and 1 FinA) from the glycerol stock, and several extinct populations in the simultaneous treatment with 2 drugs, the total number of analyzed samples was as follows: 48 FLZ, 47 FinA, 47 ANI, 48 AinF, 21 ANIFLZ and 48 YPD. The growth of the samples was visually assessed by their capacity to grow at the last drug concentration(s) after 4 × 3-day long passages in YPD medium without drugs.

Susceptibility tests

Susceptibility to flz and ani was studied in a high-throughput manner using a robot, and recording not only the endpoints but also the growth curves of all drug dilution assays over at least 18h. Susceptibility tests were performed in at least three replicates following the EUCAST E.DEF 7.3.1. protocol. Briefly, isolates were pinned on agar containing RPMI with 2% glucose buffered with MOPS (3-(N-morpholino)propanesulfonic acid) and grown at 37°C. Fresh overnight cultured strains were adjusted to 2-10 x105 CFU/mL in distilled water. Next 50 μL of broth was then added to 150 μL antifungal solution (in RPMI /w MOPS) and incubated at 37°C. OD600nm was measured every 60 - 90 min and growth was evaluated after around 18h. The range of concentrations tested was 16-0.016 μg/mL for ani 256-0.25 μg/mL for flz, following EUCAST guidelines.

Fitness and susceptibility measurements

For each sample at each drug concentration, fitness was measured as the area under the time-versus-optical density curve (hereafter referred as fAUC, calculated with the qfa package (v0.0-44 http://qfa.r-forge.r-project.org/). Minimum Inhibitory Concentration 50 (MIC50) values were calculated as the minimum concentration where the fAUC relative to the no-drug control was below 50%. If 50% of the inhibition was not met within the tested concentration range, then MIC was set to twice the maximum assayed concentration for numerical analyses in Figures 6 and S1. We also define rAUC as the area under the drug concentration-versus-relative fitness curve (AUC), normalized by the maximum AUCMAX where there is no change in fitness across the entire range of concentrations (Figure 2A). rAUC was used as a proxy for the quantitative levels of resistance for each sample. To filter out experimental artifacts, we kept the three technical replicates that were closest to the median for each sample and measure (fitness, relative fitness, MIC and rAUC).

To correct for intraspecific fitness differences, we based our fitness analysis (see Results) on the log2-ratio between the fAUC of each sample and the unevolved WT strain. This value was used as a proxy for fitness changes occurring during the experiment. Under the same reasoning, we defined strains with acquired resistance as those where the MIC was more than 2 times the WT MIC. This threshold separated our samples clearly into susceptible and resistant strains (Figures S1C and S1D). All the fitness and susceptibility measurements are in Data S1. Doubling rate per hour was inferred from the maximum slope in the time-versus-log2 (OD) data using bins of 3 time points for the analysis of EF1620_7B_ANI (see below).

Analysis of MIC and rAUC measures of antifungal drug resistance

As discussed in the main text, both MIC and rAUC measurements were correlated (Figures S1A and S1B). However, they presented several important differences that we discuss here in more detail. First of all, MIC values presented clearer increments and a bimodal distribution, making it easier to define thresholds for resistant versus susceptible samples as compared to rAUC (Figures 2B and S1D). Accordingly, we used MIC values to define resistant samples. In addition, although measurement errors are similar (Figures S1A and S1B) MIC is more consistent across independently evolved strains of the same condition (Figures 2B and S1C). However, rAUC values provided a continuous estimate of resistance, which is better suited for quantitative analyses (such as those of Figures 2E and 7C). Importantly, rAUC was not affected by the trailing effect. This effect occurs when total growth inhibition is not achieved with increasing concentration of the drug, but rather cell densities are maintained. This effect has been reported with azoles and Candida species.65, 66, 67 We observed this effect occurring in most (8/10) ANI samples with ERG3 mutations, leading to high MIC values that were in the range of FLZ samples (Figure S7C). The rAUC values, however, were not affected by the trailing effect and these strains presented flz rAUC values intermediate between flz non-resistant ANI and flz-resistant FLZ samples (Figures 2B and 7C). Conversely, there is one sample (BG2_11H_ANI) bearing an ERG3 premature termination codon and presenting a mismatch between flz MIC and rAUC. Although MIC is in the WT range, visual inspection of the flz concentration-versus-fitness curve showed a trailing effect around 50% of growth (Figure S7C), implying increased resistance. This is consistent with the observed high rAUC (Figures 7C and S7A). Taken together, these examples suggest that rAUC captures better the quantitative landscape of drug resistance.

Finally, we found another sample (EF1620_7B_ANI) where neither MIC nor rAUC captured the true nature of flz resistance. This sample shows a non-monotonic relationship between flz concentration and relative fitness (Figures S7C and S7D). This motivated us to analyze this sample under another fitness estimate, the doubling rate per hour (DR), in addition to fAUC. We found that this sample had low fitness (by both fAUC and DR) in the absence of the drug, with a small increase in the lower flz concentrations. This low level of basal fitness results in high relative fitness at low drug concentrations (as compared to other samples) (Figures S7C and S7D). This analysis suggests that this non-monotonic relationship (if present) is very weak in terms of absolute fitness. This example illustrates how MIC and rAUC values can be misleading in strains with very low basal fitness.

DNA extraction

A modified protocol from the MasterPure Yeast DNA Purification Kit was used to extract DNA. In brief, samples were grown overnight in liquid YPD at 37°C. Cells were pelleted and lysed with RNase treatment at 65°C for 15 min. After 5 min of cooling down on ice, samples were purified by the kit reagent by mixing, centrifugation and removal of the debris as described in the kit protocol. Further, samples were left at −20°C with absolute ethanol for at least 2 h after which the DNA was precipitated for 30 min at 4°C. The pellet was washed in 70% ethanol and left to dry. TE buffer was used to resuspend the DNA. The Genomic DNA Clean & Concentrator kit (Zymo Research) was used for the final purification.

Target FKS and ERG3 sequencing

All ani-exposed samples (ANI, ANIFLZ and FinA) were examined for mutations in one region of FKS1 and two regions of FKS2 encompassing echinocandin resistance mutational HSs. Three samples without mutations in the above-mentioned HSs were also inspected in the HS2 of FKS1. All the new FKS mutations are in Data S2. We used PCR primers described earlier (Table S5). ANI samples not subjected to WGS were also amplified by two PCRs with two sets of primers (Table S5) to obtain ERG3 sequences. PCRs were carried out by using Taq DNA polymerase from DongShengBio. The reaction mixture included primers of concentration of 0.4 μM, 20 μL Taq DNA polymerase, 1 μL liquid sample grown for 24-48 h in YPD and water up to a final volume of 40 μL. Optimase ProtocolWriter was used to develop conditions for each primer set.

We tested for the possible trajectories of final FKS and ERG3 mutations in the 10 ANI samples subjected to WGS and presenting ERG3 alterations to infer which might have appeared first in the evolution. We selected and analyzed single colonies from our glycerols stocks of stored populations after the 2nd passage at 0.032, 0.064, 0.128 and 0.256 ug/ml ani (beginning of the adaptation). PCRs were carried out as described above.

Petite phenotype in ani adapted mutants

10 ANI samples that underwent WGS and show changes in ERG3 gene, CBS138 WT and Saccharomyces cerevisiae petite control were inspected for presenting a petite phenotype. Samples were grown on YPD (1% yeast extract, 2% bactopeptone, 2% glucose) and YPG (1% yeast extract, 2% bactopeptone, 2% glycerol) for 24h-48h.

Whole genome sequencing

Evolved mutants: Genome sequences were obtained at the Ultra-sequencing core facility of the CRG, using Illumina HiSeq 2500 sequencing machines, and as previously described. In brief, libraries of paired-end, 125 bases-long reads were prepared. The DNA was fragmented by nebulization or in Covaris to a final size of ∼600 bp. After shearing, the ends of the DNA fragments were blunted with T4 DNA polymerase and the Klenow fragment (New England Biolabs). DNA was purified using QIAquick PCR purification kit (QIAGEN). 3′-adenylation was performed by incubation with dATP and the 3′-5′-exo-Klenow fragment (New England Biolabs). DNA was purified using MinElute spin columns (QIAGEN) and double-stranded Illumina paired-end adapters were ligated to the DNA using rapid T4 DNA ligase (New England Biolabs). After another purification step, adaptor-ligated fragments were enriched, and adapters were extended by selective amplification in an 18-cycle PCR reaction using Phusion DNA polymerase (Finnzymes). Libraries were quantified and loaded into Illumina flow-cells at concentrations of 7–20 pM. Cluster generation was performed in an Illumina cluster station. Sequence runs of 2 × 100 cycles were performed on the sequencing instrument. Base calling was performed using Illumina pipeline software. In multiplexed libraries, we used 4 bp internal indexes (5 indexed sequences). De-convolution was performed using the CASAVA software (Illumina). Sequence data of the genomes have been deposited in the Short Read Archive (SRA) database, with accession number PRJNA635652 (SRA: PRJNA635652).

The genome of the CRISPR 3H_ANI with ERG3(D122Y) sample was pooled with two genomes from divergent species (Candida albicans and Candida parapsilosis), after confirming with Crossmapper the absence of read cross-mapping in the chosen sequencing design. Sequencing libraries were made at the Functional Genomics Core Facility at the IRB and genome sequences were obtained at the sequencing core facility of the CNAG. 500-1,000 ng of genomic DNA dissolved in a final volume of 50 ul TE buffer were sheared with a Bioruptor sonicator (Diagenode) using the following settings: temperature 4-10°C; intensity: high; cycles: 3; cycle time: 5 minutes; cycle program: 30 s pulse and 30 s rest time. At the end of each sonication cycle samples were centrifuged at 4°C and the water tank was refilled with pre-cooled water. DNA fragmentation was quality controlled using the Bioanalyzer 2100 and its DNA High Sensitivity chip (Agilent) and quantified using the Qubit fluorometer and its dsDNA HS assay (Invitrogen). NGS libraries were prepared from 250 ng of fragmented DNA using the NEBNext Ultra II DNA library prep kit for Illumina (New England Biolabs). Adaptor-ligated DNA were size-selected using the provider-recommended settings to obtain an insert size distribution of 300-400 bp. After purification, libraries were amplified through five PCR cycles using the NEBNext multiple oligos for Illumina (New England Biolabs). The final libraries were quantified on Qubit and quality controlled in the Bioanalyzer. An equimolar pool was prepared with the six libraries and submitted for sequencing at the Centre Nacional d’Anàlisi Genòmica (CRG-CNAG).The libraries were sequenced on NovaSeq 6000 (Illumina) with a paired-end read length of 2x150 bp. Image analysis, base calling and quality scoring of the run were processed using the manufacturer’s software Real Time Analysis (NovaSeq 6000 RTA 3.4.4). To select the C. glabrata sequencing reads we used Burrows-Wheeler Alignment (bwa v0.7.17) mem (http://bio-bwa.sourceforge.net/bwa.shtml) to align the reads to a concatenated reference genome including the three pooled species. We took the reference genomes from the Candida Genome Database (v_s02-m07-r35 for C. glabrata and haplotype A of v_s07-m01-r110 for C. albicans) and the NCBI (sequence GCA_000182765.2 for C. parapsilosis). We next separated the reads uniquely mapping to C. glabrata with samtools (v1.9), which yielded the final whole-genome sequencing dataset.

Small variant calling and interpretation

For each library, we first performed quality control of the reads with fastqc (v0.11.8, https://www.bioinformatics.babraham.ac.uk/projects/fastqc) and trimming with trimmomatic (v0.38). The trimmed reads were aligned against the reference C. glabrata genome (the latest version by 12/03/2019, which is v_s02-m07-r35 from the Candida Genome Database (CGD: v_s02-m07-r35)) using Burrows-Wheeler Alignment (bwa v0.7.17) mem (http://bio-bwa.sourceforge.net/bwa.shtml). In addition, indexing of the genome and construction of a sequence dictionary was performed with samtools (v1.9) and picard (v2.18.26 http://broadinstitute.github.io/picard/), respectively. We next used three different algorithms (GATK Haplotype Caller (HC) (v4.1.2), freebayes (FB) (v1.3.1 https://arxiv.org/abs/1207.3907) and bcftools (BT) (v1.9, https://github.com/samtools/bcftools) to call and filter Single Nucleotide Polymorphisms (SNP) and small insertions/deletions (IN/DEL) in both haploid and diploid configurations. We defined as high-confidence (PASS) variants those with read depth above 20, with extra filters for HC and FB. For HC, we kept as PASS variants those where 1) there were less than four additional variants within 20 bases; 2) the mapping quality was above 40; 3) the confidence based on depth was above 2; 4) the phred-scaled p value was below 60; 5) the MQRankSum was above −12.5 and 6) the ReadPosRankSum was above −8. For FB, we kept as PASS variants those where 1) quality was above 1 or alternate allele observation count was above 10; 2) strand balance probability of the alternate allele was above 0; 3) number of observations in the reverse strand was above 0; and 4) number of reads placed to the right/left of the allele were above 1. We further used vcfallelicprimitives from vcflib (v1.0.0 https://github.com/vcflib/vcflib) to uniformize the called variants across the three algorithms, and the ensembl Variant Effect Predictor (v96.3) to annotate the potential functional effect of each variant in both coding and non-coding regions. In addition, we developed a tool to visualize (and better interpret) the genomic location of each variant across multiple samples using the python plotly package (v2.7). This pipeline is ready to use for any paired-end short-read sequencing library at https://github.com/Gabaldonlab/VarCall_Cglabrata_IVevolution.

We considered PASS variants to be those SNPs that passed the filtering of the three algorithms and those INDELs that passed both HC/FB filters (which were shown to have highest overlap). For each sample evolved in drug conditions, we defined variants newly-acquired during the experiment to be those that were not called in any of the corresponding WT and YPD samples. We ran this variant calling pipeline in both haploid and diploid configurations for all samples. Diploid variants may have appeared in regions that are under whole-chromosome duplications. We keep only as true “heterozygous” or “homozygous” diploid variants as those that appear to be like this by all the programs tested and within a duplicated chromosome (see below). All the new small variants are found in Data S3. In addition, Table S1 includes the variants shared between CST109 and M12 and absent in the other representatives of their clades.

Identification of large aneuploidies, segmental duplications and deletions

To detect genes affected by CNV, we calculated the read depth for each gene relative to the median read depth per gene across all nuclear chromosomes that did not have signs of large duplications (see Results) (hereafter referred to as relative coverage). The read depth was calculated using mosdepth (v0.2.6). We then defined deleted genes as those with > 50% of their length not covered by reads. To keep only gene deletions appearing during the experiment we further filtered out genes that were also lost in the corresponding WT or with a relative coverage below 0.1 in YPD-evolved sample (which may suggest a loss also in the WT or in the YPD). We manually curated the deletion list to find regions potentially deleted in a previous sample of the evolution experiment, which was the case of a small region in chromosome D (including CNE1, with a relative coverage below 0.1 in EF1620_7B_ANI) and the S. cerevisiae GPB2 ortholog (with a relative coverage below 0.1 in EF1620_7B_ANI). Importantly, these two genes were lost in a single genomic rearrangement (see below, Figure S3).

CNV was defined by calculating the log2 ratio between the relative coverage of each sample against the matching YPD (log2cov_vsYPD). Copy-number (CN) increase refers to log2cov_vsYPD above 1 and a relative coverage above 1.8, while CN decrease refers to log2cov_vsYPD below −1 and a relative coverage of the corresponding YPD above 1.8. The rationale of this filtering was to detect genes lost and under CNV during drug exposure, correcting for intrinsic biases in per-gene coverage. As noted in other studies, we found that relative coverage was correlated with the distance to the telomere (hereafter referred as “smiley-pattern”), which may be an artifact of library preparation and/or sequencing, with this effect varying across samples. We hypothesize that this is partially why most of the CNV was found in subtelomeric regions (defined here as the first and last 50 genes of a chromosome). We thus filtered out any CNV call that was not supported by equivalent genomic rearrangements (see below). In addition, chromosomes with large aneuploidies were defined as those where we consistently observe genes with increased CN and relative coverage around 2x across a region spanning at least 10% of the non-subtelomeric chromosome (Figure 5A).

Analysis of genomic rearrangements

To identify GR we implemented an algorithm that uses split-reads, discordantly aligned read-pairs and de novo assembly evidence to call genomic breakpoints and interpret the resulting GRs and CNVs. Breakpoints were called using gridss (v2.8.1) and integrated into complex structural variation with clove (v0.17). The straightforward implementation of this pipeline was challenging because of the lack of established parameters for yeast genomes, and the “simley-pattern” bias (see above) impeding the use of a single read-depth threshold for filtering deletions and tandem duplications (used by clove). We thus chose the running and filtering parameters from a simulation-based optimization implemented in the perSVade pipeline (v0.0, https://github.com/Gabaldonlab/perSVade).

GR appearing during the experiment were defined as those where none of the breakends (each of the ends of a breakpoint) matched a breakend in any of the parents (with an overlap of less 200 bp), in a way that resembles the small variant calling (see above). This is an extremely conservative approach (as most called breakends in the parents may be false positives) to ensure high confidence in our final set of variants. In addition, we defined “haploid breakends” as those with an allele frequency (AF) above 0.75 and “heterozygous breakends” as those with an AF > 0.25. We also filtered out tandem duplications, inversions and deletions where any of the breakends was not haploid, as these variants can not yield heterozygous breakends in haploid chromosomes. Note that we did not detect any such heterozygous events in the aneuploid chromosomes. Furthermore, we manually curated the results to identify errors in the summarization of breakpoints into complex rearrangements. This approach yielded one sample (P35_10E_FinA) with two reciprocal inverted interchromosomal breakpoints between close positions (less than 200 bp apart) of chromosome (Chr) G (breaking the CDS of FKS1) and ChrM. These were called as two independent unbalanced translocations, but we interpret them as an inverted balanced translocation between the two chromosomes. The coverage “smiley-pattern” was also consistent with this model.

To focus on resistance-conferring events, we examined genes with ns mutations or nearby GR (within less than 2kb) appearing recurrently (at least twice) in our experiment. These included ERG3, FKS1 and the ortholog of S. cerevisiae CNE1, mentioned in the main text (see Results). We confirmed all these rearrangements through PCR (see below). Regarding ERG3, we found one ANIFLZ sample with a deletion at the beginning of the CDS and a FinA sample with a deletion in the 5′ region (potentially spanning the promoter, and related to the loss of CAGL0F01815 g (see Results). Both of these were associated with low relative coverage (< 0.01) spanning the breakpoint, which further confirmed these deletions (Figure S3A). These are additional ERG3 mutations potentially related to ani exposure. We also found an inter-chromosomal breakpoint between ChrD and ChrL in EF1620_7B_ANI with the orientation of a deletion breakpoint. Importantly, the WT strain underwent a balanced translocation between these chromosomes (as compared to the reference genome), which means that the alteration appearing upon drug exposure was actually a deletion event (also confirmed by coverage). The deleted region included CNE1, which may be related to ani adaptation (see Results) (Figure S3B). This also constitutes an example of how the rearrangements found in each strain modulate the interpretation of breakpoints appearing during the experiment. Finally, we found two FinA samples with GR breaking the FKS1 coding region, including one deletion at the beginning of the coding sequence (with relative coverage < 0.01) and one balanced inverted translocation between ChrG and ChrM (Figure S3C). Both samples carried FKS2 mutations (potentially conferring ani resistance), suggesting that these rearrangements are complementary FKS1 alterations with a similar impact as the truncating small variants mentioned in the main text.

On another note, we attempted to infer the precise events leading to partial aneuploidies during the experiment (Figures 5A and 6). We found an unbalanced translocation explaining the partial duplication of ChrE in CBS138_9F_FLZ. Our GR-detection method predicted that the right arm of ChrE (matching the aneuploid region (Figure 5A) was duplicated and attached to ChrJ, replacing the left-end at the breakpoint. This region showed low coverage after the breakpoint (supporting the unbalanced translocation call), but not until the end of the chromosome (which would be expected from such an event). Interestingly, the deleted region was found between the unbalanced translocation breakpoint and a location with low WT coverage. We propose that this configuration is the result of a pre-existing rearrangement in the WT strain, which explains why the deleted region does not span the entire left-end of the chromosome. Accordingly, the ChrE breakend was called heterozygous, while the ChrJ was haploid (Figure S3D). Conversely, we could only find an inverted heterozygous breakpoint matching the start of the aneuploid region of ChrE in CST34_2A_AinF, which was not enough to explain the source of the duplication. Finally, we found that the (apparently) partial duplications of ChrI in the EB0911 samples are actually whole-chromosome aneuploidies. The WT EB0911 depicted balanced translocations between Chr D, I and L, generating three (mixed) chromosomes from the successive fusions. We found one of these mixed chromosomes with 2x coverage in both samples with aneuploidies (Figure S3E). Interestingly, this chromosome is much shorter than the reference ChrI, perhaps resulting in a lower fitness cost of this aneuploidy. We speculate that this is the reason why this aneuploidy is found only in this strain. Taken together, these results suggest that complex structural variation may contribute to drug resistance. They also show how breakpoint calling can explain the precise events leading to CNV and aneuploidies.

Presence of all the GRs discussed in the text was confirmed with PCR using primers specifically designed to provide amplicons only in the presence of the GR (translocations) or with a different size (deletions). Results are presented in Figure S3F. All events were positively confirmed. Primers used for each GR validation are presented in Table S5. PCRs were performed using Taq DNA polymerase from DongShengBio. The reaction mixture included primers of concentration of 0.4 μM, 15 μL Taq DNA polymerase, 1 μL liquid sample grown for 24h in YPD and water up to a final volume of 30 μL. Optimase ProtocolWriter was used to develop conditions for each primer set.

Analysis of clinical isolates’ sequencing datasets

We obtained all the variant calling files for publicly available whole genome sequences of Candida clinical isolates from the CandidaMine database (v1, https://candidamine.org, publication in progress). The MIC values for each sample were obtained by manual curation of the associated literature, when available.

In Candida glabrata, we could find these data in 126/393 clinical isolates, including resistance to fluconazole (flz 126/126), posaconazole (pos 84/126), voriconazole (vrz 91/126), isavuconazole (ivz 37/126), micafungin (mif 42/126), anidulafungin (ani 9/126) and caspofungin (cas 91/126). Some of these drugs lack established clinical resistance breakpoints, which did not allow a direct identification of resistant isolates. We thus, defined the resistance breakpoint for each drug as 2x the maximum MIC reported in a set of susceptible isolates (from Carreté et al.). Ani susceptibility was not measured for these isolates, so that we took the standard EUCAST breakpoint to define ani resistance. This data is sparse, so that we do not always know the MIC values for all drugs in a given isolate. We thus, focused our analysis on “azole” or “echinocandin” resistance instead of splitting by individual drugs. In order to achieve this, we defined an isolate to be “resistant” to a given class of drugs if it was resistant to all the measured drugs of that class. This yielded 41/126 and 19/91 isolates resistant to all tested azoles or echinocandins, respectively. We could find two samples with resistance to both classes of drugs. In Candida albicans, we could find MIC data for 187/478 clinical isolates. We could define the resistance breakpoints according to EUCAST for all tested drugs but caspofungin. We defined an isolate to be cas-resistant if the MIC was above the percentile 90. This yielded 39/186 and 9/150 isolates resistant to all tested azoles or echinocandins, respectively. We could find one sample with resistance to both classes of drugs. Given the low numbers of samples with resistance to both drugs, we conclude that the available data is insufficient to perform analysis of cross-resistance or multidrug resistance. In order to assess whether the mechanisms driving single-drug resistance in vitro are clinically relevant, we first analyzed these publicly available sequences of Candida clinical isolates. We assessed how many of the drug resistance variants described in this work were also found in these clinical isolates, which yielded little or no overlap depending on the gene (Figure S5A). We hypothesized that the underlying reason is that several mutations in the same gene can explain drug resistance (Figure 6). In order to test this we calculated the overlap between CandidaMine variants and two datasets of previously described drug resistance-mutations: the SENTRY database and a set of described PDR1 mutations from the literature.71, 72, 73 This yielded low overlaps as well, comparable to those found in our work (Figure S5A).

In addition, we inferred the expected overlap between different mutation datasets through a randomization strategy on our samples. We divided the samples carrying mutations in a given gene into two random subsets. For each subset, we calculated the number of mutations only in the subset or also found in the other subset. This process was repeated 100 times, and the results (Figure S5B) show that the overlap is comparable to the observed between datasets of different works.

We conclude that it is difficult to measure the clinical impact of the mutations described here because most of them cannot be found in the currently available isolates. However, this low overlap is expected and comparable to other datasets of well-known resistance-conferring mutations.

CRISPR-Cas9 based genetic modifications

Donor DNAs

Short fragment of ERG3 with D122Y (G364T) was ordered from Integrated DNA Technologies, Inc.. This fragment of the gene also contained additional synonymous mutations in PAM region (short NGG sequence that follows the DNA region targeted for cleavage by the CRISPR system) to bypass recutting by the Cas9 once the donor DNA is integrated, hence to improve the number of positive transformations. A large donor DNA containing ERG3 mutation (D122Y) was also amplified from 3B_ANI evolved sample by FL1_FWD and FL2_REV primers. All primers and the ordered sequence can be found in the supplementary information (Table S5). Two approaches were used to introduce ERG3 mutations. The first approach involved the transformation of a fragment containing the ERG3 alteration and assumption that the positive transformants would exhibit increased resistance to fluconazole, hence the transformation was followed by selection on agar plates containing fluconazole. Second approach involved creating a DNA construct containing ERG3 gene fused with NAT1 gene (upstream) as a selection marker. NAT1 was amplified from a vector pTS50 (a kind gift from Karl Kuchler). Two of such donors were used. One contained NAT1 fused with wild-type ERG3 (amplified from DNA extracted from wild-type Candida glabrata strain) and second contained NAT1 with ERG3 bearing the mutation (D122Y, amplified from DNA of fluconazole and anidulafungin resistant evolved mutant (3B_ANI)). The first donor was used to examine the influence of the presence of NAT1 on flz susceptibility as well to eliminate the mutation acquired during the evolution and check for the reversion of the phenotype. In this approach, ERG3 with downstream region was amplified by PCR using FLKI_ERG3 set of primers from 2 strains: one containing wild-type ERG3 and one containing the mutation (3B_ANI). Upstream ERG3 region was amplified by FLKII set of primers. NAT1 was amplified from a vector pTS50 by PCR and ‘NAT1(for DNA donor constructs))’ set of primers. All primers contain additional homologous sequences to ensure the fusion FLKI_ERG3:NAT1:FLKII. The fused fragments were gel purified and correct fusion was confirmed by PCR with internally placed primers –inside_ERG3_FWD and inside_NAT1_REV and inside_NAT1_FWD with flank_ERG3_REV.

CRISPR-based mutagenesis

CRISPR-based mutagenesis was performed using ribonucleoproteins (RNPs) and following a previously described method by Grahl et al. RNPs were created using the Alt-R CRISPR-Cas9 system bought from Integrated DNA Technologies, Inc.). The CRISPR machinery included: purified Alt-R S.p. Cas9 Nuclease V3, and guide RNA containing universal transactivating Alt-R CRISPR-Cas9 tracrRNA and target specific crRNA (Table S5).

The synthetic ERG3 fragment as well as the large donor DNA containing the ERG3 mutation were transformed into 3H_ANI sample and selected on 64ug/ml flz. The same trial of transformations was done on CBS138 WT strain but the selection of positive transformants was unsuccessful. One of the positive transformants was subjected to whole genome sequencing to infer the presence of only inserted ERG3 mutations and absence of additional protein altering mutations or CNVs (which could explain the resistance). In parallel, an alternative approach with improved selection was conducted.

ERG3 with NAT1 were transformed into wild-type Candida glabrata strains CBS138, EB0911 and its anidulafungin resistant progenies: 3H_ANI and 3B_ANI mutants. The positive transformants were selected on YPD with 200 μg/ml nourseothricin. To ensure that the DNA donors were transformed in the correct place in the genome a PCR with a REV primer that falls outside of the designed constructs and a FWD primer that falls inside NAT1 gene was performed – inside_NAT1_FWD with out_REV. Additionally, the insertion of the NAT1 was examined by amplification of longer fragment when performing a PCR with primers surrounding the place of the insertion (inside_ERG3_FWD with flank_ERG3_REV). The PCR conditions were designed with OptimaseProtocol. Presence and absence of ERG3 mutations was confirmed by Sanger sequencing.

All transformations were performed by electroporation of competent cells prepared using lithium acetate (LiAc). Overnight cultures were diluted to an optical density at 600nm (OD600) of 0.3 in 50 mL YPD and left to grow to obtain OD600 of approximately 1.6 to 2.2. Then cells were pelleted, washed once with 25ml of sterile water and resuspended in 10ml of a transformation buffer (100 mM LiAc, 10 mM Tris-HCl, 1 mM EDTA) and incubated with shaking for 1h. The cells were further incubated with shaking for 30min with 1ml of 1M dithiothreitol (DTT), washed twice with 40ml ice-cold water and once with 5ml ice-cold 1 M sorbitol before resuspension in 200 μL of ice-cold 1 M sorbitol.

CrRNAs and tracrRNA were first dissolved in RNase-free distilled water (dH2O) at 100 μM and stored at –20°C. The guide RNA was created by mixing equimolar concentrations (4 μM final) of the gene-specific crRNA and tracrRNA (to obtain a final volume of 3.6 μl per transformation) and incubating at 95°C for 5 min, followed by cooling down to room temperature. The Cas9 nuclease (60 μM stock from IDT) was diluted to 4 μM in dH2O to a volume of 3 μl per transformation. RNPs were assembled by mixing 3.6 μL of guide RNAs with 3 μl of diluted Cas9 protein, followed by incubation at room temperature for 5 min. Transformation of cells was carried out by electroporation of cell suspension containing: 40 μl of cells, 6.6 μl of RNP and 1 μg of repair constructs.

Electroporation was performed using an 0.2-cm electroporation cuvette and electroporated with a manual 1.8 pulse (Bio-Rad MicroPulser). Following the transformation, 1 mL ice-cold 1 M sorbitol was added to the cuvette. The cell suspension was then transferred to an eppendorf and the cells were gently pelleted (3 min, 3,000 rpm) before resuspension in 1 mL of YPD. Cells were recovered for 3 to 4 h at 30°C while gently shaking. After recovery, cells were pelleted and resuspended in 200 μL liquid YPD and the aliquots were spread plated onto YPD plates with 200 μg/ml nourseothricin and incubated at 37°C for 2 days.

Validation of phenotypes

Spot tests were performed to visualize changes that the transformations exert on antifungal drugs susceptibilities. Briefly, overnight cultures were set to the OD = 0.5 and serially diluted 10-fold and 10ul was spotted on YPD agar plates supplemented with antifungal drugs (Figure S2).

Fitness competition

In vitro competitive fitness was tested between ani resistant strains (containing FKS mutations) with and without ERG3 mutations and in rich medium as well as in rich medium supplemented with 0.5ug/ml of anidulafungin. To be able to distinguish the strains, we used CRISPR transformants containing NAT1 as a selection marker. Two pairs of ani resistant strains were used: 1:ERG3(D122Y) versus ERG3(WT)+NAT1 and 2:ERG3(D122Y)+NAT1 versus ERG3(WT). Two pairs were used to assure no fitness effect of the presence of NAT1. The competition test was conducted following a protocol described by Duxbury et al. Briefly, all 4 strains were grown overnight and adjusted to 6.49 × 106 cells/ml prior to mixing and subsequent two fold dilution in the growth media. The first pair of strains was mixed in two different ratios (50:50 and 75:25), while the second pair in one (50:50). Each pool of mixed strains along with the strains alone were inoculated in wells of a 96 well plate in triplicates and incubated at 37°C with shaking for 24h hours. Cells at the beginning of the experiment, after 24h growth in YPD medium as well as in YPD + 0.5 ug/ml anidulafungin were diluted and plated on YPD agar plates and YPD agar plates with 200 μg/ml nourseothricin (each at least in duplicates). The number of cells of the strains that lack NAT1 were obtained by subtracting the cells obtained from YPD+nourseothricin plate (average of the plated replicates) from the total number of cells observed on YPD plates (average of the plated replicates). Since we observed that strains containing NAT1 were growing in lower abundance on the antibiotic than on YPD plates alone, the total number of cells were accounted for this discrepancy. In vivo fitness competition between the strain containing ERG3 mutation (D122Y) and ERG3(WT)+NAT1 was also tested in Galleria mellonella. For that, the overnight grown strains were adjusted to 2.5 × 108 CFU/ml in PBS, mixed in the ratio 50:50 and 10ul of the cell suspension was injected into at least 3 larva and left 24h at 37°C. To determine the fungal burden, 3 larvas per mix were briefly washed in 70% ethanol followed by sterile water, and then placed into screw-cap tubes with 3 sterile glass beads and 1ml of PBS. The tissue was then homogenized through 3 rounds of shaking for 20 s at 4 m/s in a Fastprep-24 (MP Biomedicals). The suspensions were serially diluted and inoculated into YPD+chloramphenicol (100 μg/ml) and YPD+chloramphenicol (100 μg/ml) + nourseothricin (200 μg/ml). The real ratios of the cells at the beginning and end of the experiment were obtained as described in vitro competition experiment.

Virulence assays

The differences in virulence between strains with and without chromosomal duplications were tested in Galleria mellonella. Three strains were used: EB0911, parental WT, and its two flz evolved progenies 3B_FLZ and 3H_FLZ, where the second presents chromosomal duplications (ChrE and ChrI). Groups of 20 healthy larvae were injected with 10 μl of cell suspension, equivalent to 7.5 × 106 CFU, into the haemocoel with a Hamilton syringe through the last left pro-leg. Control set of larvae were injected with 10ul of PBS. Following infection, larvae were incubated at 37°C and survival, based on response to physical stimulation, was monitored daily for 6 days. The survival plots were created by Graphpad Prism 8.4.2.

Quantification and statistical analysis

rAUC, MIC and fitness measurements

We calculated the MIC, rAUC and fitness values for all evolved samples as explained in the STAR Methods section ‘Fitness and susceptibility measurements’. For each evolved strain and drug concentration, we measured between three to five technical replicates, and kept the three replicates that were closest to the median for each measure (fitness, relative fitness, MIC and rAUC). We used the median across these three replicates and the central estimate for several analyses (Figures 2, 6, 7, S1, and S6; Data S1). In addition, we calculated the median absolute deviation across technical replicates as a measurement of dispersion (Figures 2, 7, S1, and S6; Data S1). All these measurements were performed with python (v3.7.8) and the packages pandas (v1.1.1) and scipy.stats (v1.5.2).

Correlation analyses

We calculated the spearman correlation between fitness and drug resistance (Figure 2E), between ani resistance and the distance to the FKS hotspot (HS) in samples with no-HS mutations (see Results; Data S2), between the number of different genes with acquired mutations and fitness/susceptibility (see Results) and between the rAUC and MIC (Figure S1). All the fitness and susceptibility measurements were taken as the median across technical replicates (see the STAR Methods section Fitness and susceptibility measurements). We defined as significant correlations those with a p < 0.05. We used the python package scipy.stats (v1.5.2) to perform all these analyses.

Association between ERG3 mutations and flz resistance

We used a Fisher’s exact test to evaluate the correlation between the presence of ERG3 non synonymous mutations and flz resistance in anidulafungin-evolved samples (Figure 7D). We defined this as a significant association because of the p < 0.05. We used the python package scipy.stats (v1.5.2) to perform this analysis.

Comparing continuous distributions

We implemented a Kolmogorov-Smirnov test to assess the equality of two distributions in a non-parametric manner, thus not assuming normal distributions. We used this test to compare the relative fitness levels in each evolution condition and the YPD samples (see Results), the flz rAUC levels of anidulafungin-evolved samples with and without ERG3 mutations (Figure 7C), the intra-strain versus inter-strain differences in fitness/susceptibility (Figure S1E) and the fitness/susceptibility levels of samples with different mutations (Figure S6). We used p < 0.05 as threshold for significant differences. We used the python package scipy.stats (v1.5.2) to perform all these analyses.

Variant calling from sequencing data

We identified the genomic variants (SNPs, INDELs, CNVs and genomic rearrangements) appearing during the in vitro evolution from whole-genome sequencing data (Figures 4, 5, 6, 7, S3–S7; Table S2; Data S3) as described in the STAR Methods sections ‘Small variant calling and interpretation’, ‘Identification of large aneuploidies, segmental duplications and deletions’ and ‘Analysis of genomic rearrangements’. In addition, we used the python package scipy.stats (v1.5.2) to calculate the mean and standard deviation of the number of new small variants across replicates of the same strain and condition (Figure 4).

Competitive fitness measurements

In vitro fitness competition (Figure 7E) was performed in three replicates of each mixed ratio of the samples. Each of these mixed replicates was plated on agar plates at least twice to obtain an average number of the growing cells. Further, we calculated mean values of growing cells and standard deviation between them. In vivo fitness competition (Figure 7F) was performed in two separate competition experiments and both in three Galleria mellonella larvae. Mean number of growing cells and standard deviations were calculated from the averaged number of cells obtained from each of the homogenized larvae separately (plated on at least two agar plates). All calculations were done in Libre Office (v 6.0.7.3).

Estimating the overlap between drug resistance mutations among samples

We inferred an expected overlap between drug resistance mutations among different samples of the same condition (Figure S5B) using python (v3.7.8) and the package pandas (v1.1.1).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
Anidulafungin	CYMIT QUIMICA S.L.	Cat# 3D-FA16270-10
Fluconazole	SIGMA-ALDRICH QUIMICA S.L.	Cat# F8929-100MG
Caspofungin diacetate	SIGMA-ALDRICH QUIMICA S.L.	Cat# SML0425-5MG
Voriconazole	SIGMA-ALDRICH QUIMICA S.L.	Cat# PZ0005-5MG
Amphotericin B from Streptomyces sp.	SIGMA-ALDRICH QUIMICA S.L.	Cat# A4888-100MG
Flucytosine	SIGMA ALDRICH	Cat# PHR1659
Chloramphenicol	Merck Life Science S.L.U.	Cat# C1919-25G
Pfu Mix	DongSheng Biotech	Cat# P2022
Taq Mix, 1mlx5	DongSheng Biotech	Cat# P2012
Fluorescent Brightener 28 - Calcofluor White (1 g)	SIGMA-ALDRICH QUIMICA S.L.	Cat# F3543-1G
Congo Red	SIGMA-ALDRICH QUIMICA S.L.	Cat# C6277-25G
Hydrogen peroxide solution	SIGMA-ALDRICH QUIMICA S.L.	Cat# 16911-250ML-F
DTT, DL-DITHIOTHREITOL	Thermo Fisher Scientific	Cat# R0861
Sodium chloride, for molecular biology	PANREAC	Cat# A2942,1000
Sodium docecyl sulfate, SDS	PANREAC	Cat# A2263,0100
Methanol (Reag. Ph. Eur.) for analysis, ACS, ISO	PANREAC QUIMICA SLU	Cat# 1310911211
DMSO (Dimethyl sulfoxide), Sterile	Werfen España S.A.U.	Cat# 16712611S
MOPS	SIGMA-ALDRICH QUIMICA S.L.	Cat# M3183
RPMI-1640 (without HEPES and Sodium bicarbonate; with L-glutamine and phenol red)	SIGMA-ALDRICH QUIMICA S.L.	Cat# 51800035
DMSO (Dimethyl sulfoxide) for EUCAST	SIGMA-ALDRICH QUIMICA S.L.	Cat# W387520
Glucose monohydrate	Carl Roth GmbH + Co. KG	Cat# 6780.4
EtOH (Supelco)	MERCK	Cat# 1.00983.1011
Glycerin anhydrous/GLYCEROL 100% Molecular Biology grade	PANREAC	Cat# A2926,1000
T4 DNA polymerase	New England Biolabs	Cat# M0201L
dATP	New England Biolabs	Cat# N0440S
3′ −5′ -exo- Klenow fragment	New England Biolabs	Cat# M0212L
T4 DNA ligase	New England Biolabs	Cat# M0202L
Phusion DNA polymerase	Finnzymes	Cat# F530S
Sorbitol	SIGMA ALDRICH	Cat# S1876-500G
Tris hydrochloride	PANREAC	Cat# A3452
Lithium acetate	SIGMA ALDRICH	Cat# L4158
EDTA	SIGMA ALDRICH	Cat# E5134-500G
Critical commercial assays
MasterPure Yeast DNA Purification Kit (200 Purif.)	BIONOVA CIENTIFICA S.A.	Cat# MPY80200
Genomic DNA clean & concentrator	ZYMO RESEARCH	Cat# D4011
QIAquick PCR purification kit	QIAGEN	Cat# 50928106
MinElute PCR Purification Kit	QIAGEN	Cat# 28004
Agilent High Sensitivity DNA Kit	AGILENT	Cat# 5067-4626
NEBNext Ultra II DNA library prep kit for Illumina	New England Biolabs	Cat# E7645L
NEBNext® Multiplex Oligos for Illumina	New England Biolabs	Cat# E7335L
Qubit® dsDNA BR Assay Kit	INVITROGEN	Cat# Q32850
Qubit® dsDNA HS Assay Kit	INVITROGEN	Cat# Q32851
Deposited data
Sequence data	This study	https://www.ncbi.nlm.nih.gov/sra/PRJNA635652
Experimental models: Organisms/strains
Candida glabrata CST109	20	CST109
Candida glabrataCST 34	20	CST 34
Candida glabrata EB0911	20	EB0911
Candida glabrata CST78	20	CST78
Candida glabrata M12	20	M12
Candida glabrata EF1237	20	EF1237
Candida glabrata EF1620	20	EF1620
Candida glabrata F15	20	F15
Candida glabrata reference genome CBS138	20	CBS138
Candida glabrata P35_2	20	P35_2
Candida glabrata BG2	20	BG2
Candida glabrata SLL2 glab	This study	SLL2 glab
Recombinant DNA
vector pTS50 with NAT1 (Karl Kuchler lab)	48	pTS50
Software and algorithms
qfa package (v0.0-44), R package	49	http://qfa.r-forge.r-project.org/
Crossmapper	50	https://github.com/Gabaldonlab/crossmapper
NovaSeq 6000 RTA 3.4.4	51	https://www.illumina.com
Burrows-Wheeler Alignment (v0.7.17)	52	http://bio-bwa.sourceforge.net/bwa.shtml
samtools (v1.9)	53	http://samtools.sourceforge.net/
fastqc (v0.11.8)	N/A	https://www.bioinformatics.babraham.ac.uk/projects/fastqc
trimmomatic (v0.38)	54	http://www.usadellab.org/cms/?page=trimmomatic
picard (v2.18.26)	N/A	http://broadinstitute.github.io/picard/
GATK Haplotype Caller (v4.1.2)	55	https://github.com/broadinstitute/gatk
freebayes (v1.3.1)	N/A	https://arxiv.org/abs/1207.3907
bcftools (v1.9)	N/A	https://github.com/samtools/bcftools
vcfallelicprimitives from vcflib (v1.0.0)	56	https://github.com/vcflib/vcflib
ensembl Variant Effect Predictor (v96.3)	55,57	https://useast.ensembl.org/info/docs/tools/vep/index.html
python plotly package (v2.7)	58	https://plotly-r.com
Pipeline for small variant and CNV calling	This study	https://github.com/Gabaldonlab/VarCall_Cglabrata_IVevolution.
mosdepth (v0.2.6)	59	https://github.com/brentp/mosdepth
gridss (v2.8.1)	60	https://github.com/PapenfussLab/gridss
clove (v0.17)	61	https://www.github.com/PapenfussLab.
perSVade pipeline (v0.0)	N/A	https://github.com/Gabaldonlab/perSVade
python scipy.stats (v1.5.2)	N/A	http://www.scipy.org
Optimase Protocol Writer	N/A	http://www.mutationdiscovery.com/md/MD.com/screens/optimase/OptimaseInput.html?action=none
Libre Office (v6.0.7.3)	N/A	https://www.libreoffice.org
Graphpad Prism (v8.4.2)	N/A	https://www.graphpad.com
Oligonucleotides
Oligonucleotides used in this study—see Table S5	N/A	N/A
Other
Sandwich cover	Enzyscreeen BV	Cat# CR1296
MegaBlock 96 Well 2.2 ml Plates	Sarsted	Cat# 82.1972.002
Nunc OmniTray	Life Technologies	Cat# 242811
3mm glass beads	SIGMA-ALDRICH QUIMICA S.L.	Cat# 1040150500
Microplate, 96 well, PS, F-BOTTOM, clear, sterile, 2 PCS./BAG	Greiner Bio-One North America, Inc.	Cat# 655161
Lid, PS, High Profile (9 MM), clear, sterile, single packed	Greiner Bio-One North America, Inc.	Cat# 656161