Systematic Prediction of Antifungal Drug Synergy by Chemogenomic Screening in Saccharomyces cerevisiae

Since the earliest days of using natural remedies, combining therapies for disease treatment has been standard practice. Combination treatments exhibit synergistic effects, broadly defined as a greater-than-additive effect of two or more therapeutic agents. Clinicians often use their experience and expertise to tailor such combinations to maximize the therapeutic effect. Although understanding and predicting biophysical underpinnings of synergy have benefitted from high-throughput screening and computational studies, one challenge is how to best design and analyze the results of synergy studies, especially because the number of possible combinations to test quickly becomes unmanageable. Nevertheless, the benefits of such studies are clear—by combining multiple drugs in the treatment of infectious disease and cancer, for instance, one can lessen host toxicity and simultaneously reduce the likelihood of resistance to treatment. This study introduces a new approach to characterize drug synergy, in which we extend the widely validated chemogenomic HIP–HOP assay to drug combinations; this assay involves parallel screening of comprehensive collections of barcoded deletion mutants. We identify a class of “combination-specific sensitive strains” that introduces mechanisms for the synergies we observe and further suggest focused follow-up studies.


INTRODUCTION
Drugs and drug-like molecules are powerful molecular tools that can act by rapid and reversible inhibition of a specific protein or other biomolecule in cells. Such chemical perturbations, while similar to genetic manipulations, have several experimental advantages: they are tunable, fast-acting, often reversible, and applicable across large evolutionary distances, e.g., from yeast to human. Drugs can be easily combined to simultaneously modulate multiple proteins' activities, and in fact, the modulation of gene products by administering a combination of drugs can be vital for a successful course of treatment (Keith et al., 2005). The clinical success of chemical combination therapies has motivated our empirical study of synergistic chemical interactions. These data can then be assessed to predict how two drugs might interact in a biological system. To study the potential interactions, several mathematical models of drug synergy are available (Loewe, 1928;Bliss, 1939;Lehár et al., 2007); two widely used approaches are the Bliss model of independence (Bliss, 1939) and the Loewe additivity model (Loewe, 1928(Loewe, , 1957. Neither model is able to explain all drug synergies, and no mathematical model is suited for all observed chemical interactions, indicating the complexity of the problem. Invasive fungal infections (IFI) are lethal threats to human health, and they cause almost two million worldwide deaths annually. In 2018, the death rate among patients suffering from IFI was reported to be 28.8% (Webb et al., 2018). At present, the available therapies, particularly for invasive infections, are limited to four categories of antifungal drugs; azoles, polyenes, echinocandins, and 5-flucytosine (Perfect, 2017), and the clinical results from most IFI cases are not optimal. In addition, emerging pathogens resistant to common antifungals (Fairlamb et al., 2016) such as the pan drug-resistant yeast Candida auris have spread in health care facilities globally (Meis and Voss, 2019;Zhang et al., 2020). One potential solution to the dearth of effective treatments is to explore the antifungal efficacy of novel drug combinations, including those prescribed for diverse indications (Livengood et al., 2020). The use of drug combinations gives rise to several opportunities: (1) it has been proven, both empirically and theoretically, that drugs that are synergistic for a particular effect do not tend to show synergy for side effects (Cokol et al., 2011), (2) the dose of individual agents with serious side effects can be reduced in a combination, (3) synergistic antifungal activity increases therapy potency and reduces lengths of treatment, and (4) compared with monotherapy, it minimizes the risk for antifungal resistance (Livengood et al., 2020). Based upon The National Institutes of Health (NIH) reports on ClinicalTrials.gov, as of May 2021, there are ∼500 on-going or completed clinical trials involving antifungal drug combinations, but the success rate of such trials has been modest.
Considering the limited number of drugs available for IFI treatment (Hill et al., 2013), we sought to expand upon our strategy to use yeast as a eukaryotic model to screen any drug that inhibits the growth of, or kills yeast-Saccharomyces cerevisiae in this study. Even though such drugs may be active against the host itself, our rationale is that using these drugs could lower host-dependent side-effects because each agent in a combination is typically applied in lower doses. Screens using select combinations on Candida spp. and Saccharomyces have been performed, but in this we work cast a broader, unbiased net (Hill et al., 2013;Shekhar-Guturja et al., 2016). Specifically, in this study, we selected 11 compounds based upon their well-characterized targets in yeast, and among all possible combinations, drug pairs that empirically showed synergy were used in HIP-HOP assays-a validated genomewide screen based on HaploInsufficiency Profiling (HIP) and HOmozygous Profiling (HOP) to quantify the relative abundance of uniquely tagged yeast deletion strains. Briefly, in the rapid and cost-effective HIP assay, a complete collection of heterozygous deletion strains is pooled, grown in the presence of the compound, and sampled as a function of time. Molecular barcodes incorporated into each strain allow parallel analysis and relative strain abundance to be quantitatively assessed either by hybridization to oligonucleotide arrays, or more recently, by Next-Generation Sequencing Technologies. The result is a list of genes ranked in order of their importance for growth and survival, a quantitative metric termed "fitness." Strains most sensitive to drugs often carry deletions in genes that encode the drug target. Once the primary mechanism has been identified and confirmed in secondary genetic and/or biochemical assays, further pathway-specific genes that act to buffer the drug target pathway can be uncovered using our HOP assay. This assay identifies the drug effects on the non-essential fraction of the genome and reveals the genes important for buffering the drug target pathway. As described in this study, we modified the HIP-HOP assay to study drug combinations. Using this approach, we were able to: (1) identify numerous synergistic combinations, (2) quantify this synergy and identify combination-specific sensitive strains on a genome-wide scale. This study is intended to generate a new strategy to predict drug synergy by using comprehensive genome-wide screens. None of the drug combinations identified here were actually tested in a fungal pathogen. However, we are hoping that this study will inspire research on relevant fungal pathogens.

Pair-Wise Screening of Drug Combinations
To identify synergistic combinations, one should determine the effect of both individual agents and drug combinations on yeast growth rate. To accomplish this, we screened the drugs in a checkerboard matrix, in which, along each axis, one of the drugs is added at progressively higher doses. Drug concentrations were selected based on inhibitory concentration (IC), which was determined prior by prescreening the drugs' effects on the wildtype cell growth; concentrations are such that there is an IC0 (no drug), IC2, IC5, IC10, IC20, and IC50 for each drug in the matrix. Hence, each drug pair was screened in a 6-by-6 dose-response curve at the IC-values listed above. All possible pairs of the 11 drugs were screened, resulting in 55 combinations, including self-by-self drug combination. The diploid yeast strain BY4743 (derived from S288C) was used to screen these combinations and was grown in a 96-well microtiter plate at 30 • C for 24 h in a TECAN optical density reader. The optical density (OD 600 ) was measured at 15-min intervals. A diagram depicts how the combination screening was performed (Figure 1).

Determination of Synergistic Combinations Based on Growth Curve Analysis
Following growth, data for all growth curves were extracted using AUDIT software (Coutin et al., 2020) as described. First, the curves were smoothed, and the area under the curve was calculated. The area under the curve was then compared to the area of no drug control (AREA drug /AREA control ) to create an inhibition ratio. We then used the Bliss multiplicative model 3 to calculate epsilon for each dose matrix, ε = Drug AB ratio -(Drug A ratio × Drug B ratio ). Specifically, we considered "drug epsilon" to be the difference between the actual combined growth and the "expected" from the multiplication of the two single agents. For example, if drug A grows at 90% compared to no drug and drug B grows at 80% compared to no drug, the expected defect FIGURE 1 | Combination screening and transformation of growth data into a quantitative metric. (A) Illustrates the combination of compounds cerivastatin and fenpropimorph as an example. A 6-by-6 dose matrix is screened with an increase of each drug along the x and y axes. The blue triangles represent increasing concentration of drugs. In (B) each square in the 6-by-6 matrix is represented by a growth curve for each drug combination, which is optical density (OD 600 ) vs. time. Panel (C) shows an example of greater than additive inhibition due to synergy. The blue growth curves represent the DMSO (no drug) control. Green curves show growth in the presence of compounds A and B. The combination of A and B is shown in the bottom right-hand cell in which the black line indicates the expected growth rate based on the multiplicative model, while the red curve is the actual cell growth in the drug combination. In the heatmap (D), the color of each square reports the epsilon value generated using the multiplicative model, where black represents no interaction, yellow represents synergy, and blue represents antagonism. In this heatmap, the 6-by-6 dose response matrix of cerivastatin vs. fenpropimorph provides an example of synergy, because epsilon is negative (yellow). W AB represents the relative growth of yeast in the presence of both drug A (cerivastatin) and drug B (fenpropimorph) when compared to a no-drug control (DMSO), W A represents the relative growth of yeast in the presence of drug A when compared to a no drug control. The color-coded scale bar from yellow (synergy) to blue (antagonism) covers the spectrum where W ab < W a × W b , to no interaction W ab = W a × W b to antagonism W ab > W a × W b .
would be 90 × 80% (e.g., 72%). If the actual combination grows at 50% compared to no drug then epsilon would be 50-72% = −23% When epsilon is zero, then no interaction is observed; when epsilon is negative, there is a synergy, and positive epsilon denotes antagonism. An antagonistic interaction indicates that one of the drugs buffers the effect of the second agent- Figure 1 illustrates how growth data was transformed into a quantitative trait to determine epsilon.

Drug-Gene Interaction Screening Using Isogenic Cultures
To determine whether a deletion mutant was hypersensitive to the drug, we had to know the growth rate of both the mutant and wild-type strains with and without drug. We used heterozygous deletion mutants of the known drug-targets (listed in Table 1).
The yeast strains BY4743 and corresponding heterozygous mutants were grown, as isogenic cultures, in 96-well microtiter plates, at 30 • C for 24 h in a TECAN optical density reader. The optical density (OD 600 ) was measured at 15-min intervals. Here, the growth metric average generation time (aka AvgG) was used to assess the fitness of wild-type and mutant strains with and without drug; this metric is comprehensively described in the protocol written by our lab (Proctor et al., 2011). The generation time is calculated from the initial inflection point at the end of log phase until the second inflection point as the culture reaches stationary phase. We normalized each strain's fitness to the wild-type and subtracted any single mutant fitness that was contributed by any particular mutant, i.e., we normalized the various heterozygous mutants' growth to wildtype to take into account any fitness defect that was caused by haploinsufficiency.  (Turi and Loper, 1992;Truan et al., 1994) HIP rapamycin Tor2 (Heitman et al., 1991;Sabatini et al., 1994) MSP/HIP cantharidin Glc7 (Li and Casida, 1992;Honkanen, 1993) HIP fenpropimorph Erg24 (Marcireau et al., 1990;H. Lai et al., 1994) HIP latrunculin A Act1 (Spector et al., 1983;Yarmola et al., 2000) Resistance mapping benomyl Tub1 (Davidse and Flach, 1977;Sheir-Neiss et al., 1978) HIP and mutant mapping sodium fluoride Ipp1 (Yan et al., 2008) HIP hydroxyurea Rnr1 (Elledge and Davis, 1990;Tai et al., 2017) HIP and resistance mapping * Multicopy suppression profiling (target overexpression).

Predicting Synergistic Combinations via Chemogenomic Interactions
Following up on how drug-drug interactions predict druggene interactions -to predict synergy using chemogenomic data, we examined 18 datasets (see Table 2) and assessed if the known drug-targets listed in Table 1 were sensitive in any of the treatments based on the log 2 ratio of control over treatment.
We then identified combinations available in our laboratory for testing; 25 combinations were selected based on the availability of compounds present in the Giaever/Nislow laboratory as well as Boone lab, at the University of Toronto. To determine if this method can successfully predict synergistic combinations, the chances of observing synergy between randomly paired compounds need to be known.

Determination of Background Synergy Rate and Experimental Validation of Predicted Combinations
To define enrichment for synergistic combinations, the chances of observing synergy between randomly paired compound combinations must be known. To address this, 105 combinations in a 4-by-4 dosage matrix were screened. We used a smaller matrix in this experiment to maximize the number of combinations that could be screened in a short time. As a result, six combinations can be screened per 96-well plate per TECAN plate reader, instead of two combinations when screened in a 6by-6 matrix. The drug concentrations were such that there was an IC0 (no drug), IC10, IC20, and IC50 for each drug in the matrix ( Table 2). The yeast strain, BY4743, was used to screen these combinations and was grown in a 96-well microtiter plate at 30 • C for 24 h in a TECAN optical density reader. The optical density (OD 600 ) was measured at 15-min intervals.

Pooled Competitive Growth Assays
Two deletion pools, a homozygous deletion pool of 5054 strains representing non-essential genes and a heterozygous pool of 1,194 strains representing genes essential for viability, were thawed and diluted in YPD to an OD 600 of 0.0625; 700 µl cultures were then grown at 30 • C with a chemical inhibitor(s) applied at a dose that produced 10-20% growth inhibition of wild-type. An automated liquid handler robot was used to maintain logarithmic growth of pools by collecting 0.7 OD 600 s of heterozygous pool following 20 generations of growth, and 1.4 OD 600 s of homozygous pool following five generations of growth, for further processing as described below.

Assessing Fitness of Barcoded Yeast Strains by Barcode Microarray
Except where indicated, pooled assays were performed as described in the protocol by Pierce et al. (2007). Genomic DNA was isolated from cells and barcodes, amplified, and hybridized to barcode microarrays, where each barcode deletion mutant is represented by 10 hybridization signals (the uptag and downtag for each strain are each represented on the array five times; Pierce et al., 2007). Array measurements were quantile normalized such that all tags hybridized with the sample pool had similar distributions. Following normalization, we applied a correction factor to the array data to correct for feature saturation (Pierce et al., 2007) and determined the fitness of each barcoded deletion strain using the average of both barcodes. A Z-Score was calculated based on the average barcodes signal intensity against a control probe sets distribution. Positive fitness defect scores signify a decrease in strain abundance during drug treatment.

Haploinsufficiency Profiling and HOmozygous Profiling of Synergistic Combinations
A key parameter in performing genome-wide screens in yeast is to determine an appropriate screening dose. This value has been empirically determined to be 10-30% of inhibition of wild type growth (Lee et al., 2014). In practice, when performing synergy screens with two agents, one must eliminate any effects due to the action of a single agent alone. Therefore, we screened each single agent at its IC20, as well as at the dose that was used to generate a combined IC20. We therefore needed data from both the combination and individual agents. For each combination genome-wide assays five screens were performed. Specifically, Agent A at its IC10-30, Agent A the dose used when combined (usually an IC2), along with Agent B at its IC10-30, Agent B at the dose used when combined (usually an IC2). Accordingly, each genome-wide synergy assay comprises five separate screens: (1) the combined screen (A+B), (2, 3) each single agent at an IC10-30 (A and B), (4, 5) each single agent at the doses used in A+B.

Analysis of Combination Profiles
Defining the Sensitivity Score For each drug treatment, a Z-Score based on the averaged Up and DOWN barcode signal intensity was calculated. Using a Z-Score of >2, we defined a list of sensitive strains in each treatment. Filtering out strains that arose from the single drug treatments, we were able to identify unique, combination-specific sensitive strains. To aid in the analysis of drug combinations, we defined factor ε that is the sensitivity of genes in the combination minus the sum of sensitivity in the single agents-ε = Z-Score AB -(Z-Score A × Z-Score B ).

Clustering of Combination HIP-HOP Profiles Sensitivity Scores
We took raw intensity values from the barcode microarrays and normalized the logged raw intensities using a method called Supervised Normalization of Microarrays (SNM). Supervised Normalization of Microarrays (Mecham et al., 2010) was supplied with batch definitions-each batch as the arrays that have the same chip date. We also supplied this method with array descriptors corresponding to the chemical treatments (i.e., which compounds were used); this is meant to preserve biologically relevant signal. After SNM, we selected either the uptag or downtag for each strain based on the lowest variation coefficient to avoid noisy tags. The logged intensities were then used to compute robust Z-Scores. We used the median and median absolute deviation to calculate the Z-Score, and then clustered strains and chemical treatments separately. The similarity between strains/chemical treatments was based on the Pearson correlation of Z-Scores.

Examination of Gene-Ontology Terms
Gene-Ontology has the aim to standardize terms for describing gene products. This vocabulary defines a set of cell terms for which a gene can be annotated to. These annotations cover a vast range from location within the cell to specific cellular functions such as nucleotide excision repair. In this study, we used Gene-Ontology terms with more than five genes and <200 genes. To determine enrichment, we used the sensitivity score. Following ranking each gene sensitive in a specific combination, we used Gene Set Enrichment Analysis (GSEA; Subramanian et al., 2005) to determine enrichment in each category.

Visualization of HIP-HOP Screens
To facilitate the visualization of the single agent and combination screens, we provide a custom shiny app to (1) upload excel or text files, (2) visualize each screen as a scatterplot (HIP and HOP plots side-by-side), (3) rapidly identify combination specific strains (red), significantly sensitive strains (green), and non-significantly sensitive strains (violet), and (4) the significant strains are detailed below the plots with hyperlinks. Figure 4 provides a static view of the output of this app which can be accessed here: https://ggshiny.shinyapps.io/GOappCN/ Additional features allow for customization, including selecting axis limits, label font sized, thresholding for significance, and sizing of data points.

Synergy Screens
Using our database of 3,200 drugs and drug-like molecules (Lee et al., 2014) we selected 11 compounds ( Table 1) and screened all possible combinations of these drugs in a 6by-6 dosage matrix using growth curve analysis (Figure 1). The growth data were analyzed, using a computer program (AUDIT; Coutin et al., 2020) that converts raw absorbance values into growth curves. Next, to determine potential synergy, corresponding growth curve data were examined to produce the average epsilon score for the drug combination 6-by-6 matrix (defined as AvgS) from which we generated heatmaps ( Figure 1D; drug screens, heatmaps, and MATLAB code are available in Supplementary Data Files 1 and 2). If two individual drugs act independently, their effects are expected to combine multiplicatively. In other words, if a drug affecting gene x causes a fitness effect W x , and a drug affecting gene y causes a fitness effect W y , then the total effect of the drug combination (W xy ) is predicted to be W x × W y . For our purposes, we measured the deviation epsilon (ε) from this expectation (where ε xy = W xy -W x × W y ) (Onge et al., 2007;Díaz-Mejía et al., 2018). Using this score and a threshold of AvgS < −0.05 to score a combination as synergistic, 33% of all combinations showed synergy. "Large-scale experimental drug synergy screens have found that synergistic drug pairs are rare (4-10%) (Cokol et al., 2011)." Hence, given this unexpectedly high level of synergy, we applied additional filters such as additivity (Loewe, 1928(Loewe, , 1953(Loewe, , 1957, highest single agent, and a potentiation model . Finally, we identified 10 combinations that deviated from expectation in all 4 models and were therefore classified as synergistic. Interestingly, fenpropimorph vs. miconazole, fenpropimorph vs. cerivastatin, and miconazole vs. cerivastatin all possessed a strong synergy when combined. These are the only compound pairs that target the same pathway-consistent with the idea that drugs targeting the same essential pathway can be an effective means to produce synergistic combinations.

Using Synergistic Drug Combinations to Predict Drug-Gene Interactions
Having identified a high-confidence set of synergistic drug-drug interactions, we tested if these interactions could be recapitulated by combining the relevant drug-gene interactions. Our rationale was that if two drugs were synergistic, a loss-of-function mutation (as exemplified in the heterozygous state) in one of the known drug-targets would confer hypersensitivity to the second compound. In other words, if Drug A targets protein A and Drug B targets protein B to produce synergy, then one would predict that Drug A when combined with a loss-of-function mutant in B, should phenocopy the drug combination (it is worth to mention that loss-of-function screens of heterozygotes assume they are acting as recessive alleles because they do not manifest a growth defect in the absence drug). In our experience, all the heterozygote diploids behave as loss-of-function alleles with respect to their response to drug. Because each deletion allele is grown competitively along with ∼6,000 strains that are wildtype at the locus they behave as recessive alleles. To empirically test this prediction, we selected heterozygous deletion mutants of the known drug targets and challenged them with each drug listed in Table 1, and the results were compared to data derived from the combination screens.
Eleven heterozygous deletion mutants-each deleted for one of the known drug-targets-along with a wild type control strain were profiled in each drug, resulting in 121 druggene interaction tests (i.e., 11 drugs against 11 heterozygous deletion mutants). The act1 heterozygote displayed a significant fitness defect without drug treatment and was eliminated from further analysis. The remaining 110 drug-gene interactions were examined for drug sensitivity using a cut-off of ≥10% fitness defect (i.e., an inhibitory concentration of 10 or IC10) and identified 76 negative drug-gene interactions (bold numbers in Table 2) among them, 10 are the expected HIP-drug or genedrug interactions (i.e., a drug-target mutant being sensitive to a drug known to target the product of that locus; red cells in Table 2), while the remaining 66 negative drug-gene interactions were novel. Among the 66 negative drug-gene interactions, 17 of the 18 predicted interactions (yellow cells in Table 2) showed a fitness defect ≥10%, giving a significant enrichment (p-value = 0.003) of drug-gene interactions, which are predicted by synergistic combinations. Using a more stringent cut-off of greater than a 30% fitness defect, we found 23 negative drug-gene interactions (underlined numbers in Table 2), 8 of which were from the predicted drug-gene interactions (p-value = 0.0002); Table 2. Situations where two loss-of-functionscaused by either mutation or drug inhibition-leads to a more deleterious effect than the fitness reduction expected from the combined loss of individual genes are referred to as negative or aggravating interactions (for instance, synthetic sickness, or synthetic lethality; Beltrao et al., 2010).

Random Combinations of Compounds Are Rarely Synergistic
To put our drug-drug synergy observations in context, we sought to determine the chances of observing synergy when two randomly selected compounds were combined. In other words, establishing the likelihood observing synergistic effects of any two compounds would then allow one to calculate any enrichment over random chance. Accordingly, we tested all pairwise combinations of 15 random compounds in a 4-by-4 dosage matrix (Table 3). Our criteria for selecting such "random compounds" included (i) they were bioactive in yeast, and (ii) their HIP-HOP profiles showed a similar number of sensitive strains when compared to the compounds in Table 2. We further evaluated these compounds by mapping them onto the synthetic genetic array (SGA) network of gene-gene interactions (Costanzo et al., 2010). Using the multiplicative synergy model, 17% of these random combinations were synergistic (ε < −0.20) which dropped to 9.5% of combinations when the overrepresented compounds that affect the cell wall and secretion were excluded. This value is similar to previously reported combination screening studies, which report ∼10% baseline synergy in any combination (Yeh et al., 2006;Lehár et al., 2007;Farha and Brown, 2010). Figure 2 represents an example of synergy prediction screens in a 4-by-4 matrix (all the SGA heatmaps are available in Supplementary Data File 3).

Predicting Synergy
"Chemical space-a term used to encompass all possible small (>500 atoms) organic molecules, including those in biological systems-is vast" (Odling-smee and Dobson, 2004). Furthermore, the current purchasable, and readily screenable compounds comprise approximately 8 million unique compounds (Chuprina et al., 2010). Screening each of these molecules as single agents is quite daunting, while screening all pairwise combinations is impossible. The number of potential combinations is (N 2 -N)/2, with the number of starting compounds being (N). To expedite the identification of synergistic combinations, we test if we could uncover synergistic combinations from combination chemogenomic data.
We reasoned that if a drug induces a fitness defect in a particular gene-deletion mutant, but does not directly inhibit that gene product, then this drug might be synergistic when combined   with a second compound that does inhibit that gene product. To survey the possible drugs and mutants that satisfy these criteria, we first used our lab database of several thousand chemogenomic assays (Lee et al., 2014) to define when a heterozygous deletion mutant of a known drug-target is sensitive. We then sought to uncover synergistic interactions so the drug can be paired with a second drug that inhibits the known drug-target (Figure 3).
This principle is illustrated with miconazole, an antifungal that targets the enzyme Erg11. Because the hmg1 deletion strain is also sensitive to miconazole, we hypothesized that the combination of miconazole and an HMG1 inhibitor (e.g., cerivastatin), would be synergistic. To test this hypothesis, we first examined the dataset in reference (Hillenmeyer et al., 2008) and all single-agent screens performed in our lab to ask if FIGURE 3 | Diagram of synergy prediction method. The HIP-HOP profile for drug A, in which red circles are essential genes and blue are nonessential genes. Genes are listed on the x-axis and sensitivity to drug A is on the y-axis. The most sensitive strain in this example is the target of drug A. However, another essential gene also displays sensitivity to drug A. This gene is the known drug-target of drug B. Using the synergy prediction method, drug A and drug B would be predicted to be synergistic. By testing them in a drug dosage matrix, one can determine if this assumption is true.
any of the drug targets in Table 1 exhibited a fitness defect (Supplementary Data File 4). From this survey, 26 predicted synergistic combinations were selected and empirically tested to determine if their synergy rate was greater than the background synergy rate of 17%. We found ∼50% of the tested pairs were synergistic (ε < −0.20): a 3.1-fold significant enrichment over random pairs (p-value < 0.01). This approach is conceptually distinct from another synergy prediction method introduced by Jansen et al. (2009) (a bioinformatics-based approach predicting antifungal synergy using chemogenomic profiles to identify compound profiles that have a statistically significant degree of similarity to a fluconazole profile). In contrast, our approach is empirical, asking if a drug, selected based on its drug-gene interactions in HIP-HOP, can induce synergy.

HIP-HOP Combination Profiles
We next used a variation of the HIP-HOP assay to screen some of the newly identified synergistic drug combinations, testing compound combinations genome-wide. Based on these tests, all 14 confirmed synergistic combinations and 12 non-synergistic combinations, and each single agent was used for chemogenomic screening (Supplementary Data File 5). We then used the chemogenomics approach described by Lee et al. (2014) to identify any significantly sensitive strain in all screens. To further scrutinize the drug combinations, we defined the epsilon value (ε) as a fraction of uniquely sensitive genes, in both synergistic and non-synergistic pairs. We defined uniquely sensitive genes as those that are significant only when both compounds are tested together at a fitness defect score of 2.0 or greater. Two examples are shown in Figure 4, and the entire dataset of single agent and combination HIP-HOP profiles which is visualizable in our shiny app is in Supplementary Data File 4 and 5.
Although the primary goal of these genome-wide combinations is to serve as a resource for focused tests of individual combination-specific genes, several high-level observations are noteworthy: (i) combinations vary greatly in the number of specifically sensitive genes, (ii) in some cases the combination-specific strains appear to be subject to potentiation by one of the two agents (i.e., these strains can be detected at higher doses of the single agents , and (iii) the combination-specific genes identified are consistent with known mechanisms of actions of one or both of the drugs used in the combination.
We examined each combination screen, both synergistic and non-synergistic, and examined the biochemical pathways enriched in each pathway. Examining the Gene Ontology (GO) enrichments via the synergy score (Figure 5), we found that each combination provides a unique signature. For instance, the miconazole-cerivastatin combination screen was enriched for gene deletion strains involved in cell wall, cytokinesis, vesicle-related processes, and sterol biosynthesis. In contrast, the miconazole-hydroxyurea screen is enriched for vesicle-related processes and cytokinesis but not for cell wall-related or sterol biosynthesis processes. These combination-enriched GO terms can identify which cellular processes are providing resistance to the combination and could help to understand the mechanism of synergy on a combination-specific basis.

DISCUSSION
In this study, we used genome-wide chemogenomic profiling to select drug combinations for synergy testing then confirmed our predictions using combined chemogenomic assays. An interesting observation from the drug combination data is that the three inhibitors affecting the ergosterol pathway were highly synergistic when applied in combination, suggesting that compounds that inhibit different points within a pathway are more likely to be synergistic, consistent with  and the observations by (Cokol et al., 2011) that found similar compounds to be "promiscuously" synergistic. We further demonstrated that, among the synergistic combinations, 78% of the combinations tested (Table 1) were the result of combining an ergosterol inhibitor with a second agent. This indicates that the ergosterol inhibitors are highly synergistic with other agents, which is likely due to their effect on the yeast cell membrane, thereby allowing compounds more effective entry (Farha and Brown, 2010).
Using a modified genome-wide assay we demonstrated that synergistic combinations result in uniquely sensitive strains that are specific to the combinations and are not observed in either FIGURE 4 | HIP-HOP screening of select drug combinations. For each combination screen, a drug combination that inhibited wild-type yeast growth by 20% was selected and screened alongside each single agent. To identify combination-specific strains, we required that the fitness defect in the combination be 2.0 or greater and that in each single agent, the fitness defect for that strain did not exceed 2.0. For clarity, the heterozygous essential (HIP scatterplots) and homozygous non-essential (HOP scatterplots) data are plotted separately. Significantly sensitive strains are highlighted in green, and combination-specific strains are depicted in red, while violet shows non-significantly sensitive strains. The fitness defect scores are shown on the y-axes. In the case of the fenpropimorph vs. rapamycin combination (top plots), only a single strain, the essential gene DML1, was identified as a combination-specific strain. This gene product has been implicated in diverse aspects of mitochondrial function. For the hydroxyurea vs. cerivastatin combination (bottom plots), a larger number of combination-specific strains are apparent. Among these are essential genes involved in sphingolipid biosynthesis (LCB2), mitochondrial metabolism (SDH3, JAC1) as well as cell cycle checkpoints, and protein degradation at the metaphase anaphase transition (LCD1, CDC23, and CBF2). Non-essential strains specific to this combination include those involved in response to diverse stresses (STE3, CGR1) and targeted protein degradation (PEP5, KEX2).
of the single agents. Because we used stringent cut-offs, the difference we found between synergistic and non-synergistic combinations likely represents a minimum level of enrichment. We also found that each combination has its own pathway and GO enrichments (Figure 5).
During the course of this work, we confirmed that these drugdrug interactions (derived from drug combination treatment) can be recapitulated using drug-gene interactions by directly assaying loss-of-function (heterozygous deletion) mutants for a drug's known target with a drug that inhibits a synergistic target. We further found that drug-gene interactions derived from synergistic drug-drug interactions were enriched for negative interactions. To extrapolate these observations, we analyzed our single-agent chemogenomic screening data to predict combinations that might exhibit synergy. Given that we observed the baseline level of synergy between 10 and 17%, between 83 and 90% of any random combination should not be synergistic. Our approach reduces the screening required by at least 3.1-fold. This experimental approach involves: (1) using chemogenomic data to identify drugs able to make known drugtargets haploinsufficient, (2) pairing the strain with the expected drug, and (3) screening a dose matrix for synergy. In a pilot of 26 combinations, we identified 14 synergistic combinations. This method is easily adaptable to include new drug-targets, as we limited our search to 11 well-characterized drug-targets in Table 1 and only examined one dataset (Hillenmeyer et al., 2008).
Synergistic effects (between either genes or drugs) have received renewed attention, especially in light of increasingly sophisticated computational approaches and the precise genome engineering possible with CRISPR-based technology.
FIGURE 5 | Combination-specific Gene-Ontology (GO) enrichment HIP-HOP data. Examination of Gene-Ontology enrichment using ε score. In the clustered heat map generated from the ε scores, on the x-axis we have each drug combination listed on the bottom; at the top on the x-axis is a purple box which denotes if a combination is synergistic or not. GO terms are denoted on the y-axis. Red shows a particular GO term is highly enriched in the combination, gray denotes no enrichment, and blue shows significant enrichment among genes with low ε scores.
For example, Cokol et al. (2018) developed a computational framework called MAGENTA to investigate the impact of microenvironment on antibiotic combinations, stating that it enables tailoring antibiotic therapies based on the pathogen microenvironment. To predict synergistic or antagonistic interactions on various microenvironments, MAGENTA leverages chemogenomic profiles of both single drugs and metabolic perturbation. They reported several synergistic combinations against Escherichia coli and A. baumannii, and predicted bactericidal drug-combinations' effectiveness when grown in glycerol media and classified genes in the glycolysis and glyoxylate pathways as top predictors of synergy and antagonism, respectively (Cokol et al., 2018).
In 2016, Wong et al. leveraged combinatorial genetics en masse (CombiGEM) to systematically study gene and drug combinations modulating biological phenotypes (Wong et al., 2016). Combi-GEM allows for the rapid construction of barcoded, combinatorial genetic libraries that can be quantified with high-throughput sequencing. They applied CombiGEM-CRISPR to generate a library of 23,409 barcoded dual guide-RNA (gRNA) combinations, performing a high-throughput pooled screen to find gene pairs that combine to inhibit ovarian cancer cell growth. In the same study, small-molecule drug pairs were tested against the pairwise synthetic lethal hits, revealing that they exert synergistic antiproliferative effects against ovarian cancer cells (Wong et al., 2016).
Combining chemogenomics and genetic interactions, Weinstein et al. (2018) studied antifungal combinations applied to two yeast species, C. albicans and S. cerevisiae. This study showed that, both synergistic and antagonistic combinations increase the cell-type selectivity of growth-inhibiting drugs. The authors speculate that drug interactions might shift selectivity in comparison to single-drug effects in mixed microbial communities. Indeed, few drugs or drug combinations should be expected to encounter the idealized conditions in laboratory experiments-the variations observed by Weinstein et al. (2018) can change the selectivity of compounds, i.e., inverting, diminishing, or enhancing therapeutic windows.
In a recent CRISPR/Cas9 screen Huang et al. sought to identify genes whose depletion causes synthetic lethality with the broad-acting but not particularly potent Aurora kinase inhibitor VX-680 (Huang et al., 2020). They reported that HCT116 cells showed hypersensitivity to VX-680 when Haspin-a serine/threonine-protein kinase encoded by the GSG2 gene-was either depleted by CRISPR knockout or with Haspin inhibitors, confirming the synergistic effect between VX-680 and Haspin depletion or inhibition (Huang et al., 2020). Recently, Zhou et al. reported a CRISPR-based, multi-gene, knockout screening system for assembly of barcoded, high-order combinatorial guide RNA libraries, en masse. Although combination therapies promise to improve treatment efficiency of various diseases, only a few effective drug combinations-especially those employing three or more drugs (see Table S1 in reference Zhou et al., 2020)-have been introduced so far. Zhou et al. used this approach to systematically identify both pairwise and three-agent synergistic therapeutic target combinations. Their study claimed to uncover double-and triple-combinations that suppressed cancer cell growth and afforded protection against Parkinson's disease-associated toxicity (Zhou et al., 2020).

CONCLUSION
In this work, we introduce a strategy to use comprehensive genome-wide screens to first predict compounds that might be synergistic and then test novel combinations empirically. This approach should be extensible to other models and allow for a rational approach to selecting effective drug combinations. Though none of the drug combinations identified here were tested in a fungal pathogen, we hope this study sheds some light on this research field and would inspire other scholars work on relevant fungal pathogens.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
CN and AS designed the experiments. AS performed the experiments. CN, HG, and AS wrote the paper. GG supervised the data analysis. All authors analyzed the data, contributed to the article, and approved the submitted version.