The Influence of Heteroresistance on Minimum Inhibitory Concentration, Investigated Using Weak-Acid Stress in Food Spoilage Yeasts

ABSTRACT Populations of microbial cells may resist environmental stress by maintaining a high population-median resistance (IC50) or, potentially, a high variability in resistance between individual cells (heteroresistance); where heteroresistance would allow certain cells to resist high stress, provided the population was sufficiently large to include resistant cells. This study sets out to test the hypothesis that both IC50 and heteroresistance may contribute to conventional minimal inhibitory concentration (MIC) determinations, using the example of spoilage-yeast resistance to the preservative sorbic acid. Across a panel of 26 diverse yeast species, both heteroresistance and particularly IC50 were positively correlated with predicted MIC. A focused panel of 29 different isolates of a particular spoilage yeast was also examined (isolates previously recorded as Zygosaccharomyces bailii, but genome resequencing revealing that several were in fact hybrid species, Z. parabailii and Z. pseudobailii). Applying a novel high-throughput assay for heteroresistance, it was found that IC50 but not heteroresistance was positively correlated with predicted MIC when considered across all isolates of this panel, but the heteroresistance-MIC interaction differed for the individual Zygosaccharomyces subspecies. Z. pseudobailii exhibited higher heteroresistance than Z. parabailii whereas the reverse was true for IC50, suggesting possible alternative strategies for achieving high MIC between subspecies. This work highlights the limitations of conventional MIC measurements due to the effect of heteroresistance in certain organisms, as the measured resistance can vary markedly with population (inoculum) size. IMPORTANCE Food spoilage by fungi is a leading cause of food waste, with specialized food spoilage yeasts capable of growth at preservative concentrations above the legal limit, in part due to heteroresistance allowing small subpopulations of cells to exhibit extreme preservative resistance. Whereas heteroresistance has been characterized in numerous ecological contexts, measuring this phenotype systematically and assessing its importance are not encompassed by conventional assay methods. The development here of a high-throughput method for measuring heteroresistance, amenable to automation, addresses this issue and has enabled characterization of the contribution that heteroresistance may make to conventional MIC measurements. We used the example of sorbic acid heteroresistance in spoilage yeasts like Zygosaccharomyces spp., but the approach is relevant to other fungi and other inhibitors, including antifungals. The work shows how median resistance, heteroresistance, and inoculum size should all be considered when selecting appropriate inhibitor doses in real-world antimicrobial applications such as food preservation.

R esistance of microorganisms to inhibitory agents is a growing problem in therapeutics (1,2) and food supply chains (3). The resistance of a microorganism to a stressor is typically described using a conventional MIC measurement, which reflects the concentration of a stressor required to prevent growth at a particular fixed inoculum size (Box 1). Other measures of stress resistance include IC 50 , a population-median measure describing the concentration required to inhibit either 50% of cells in a population (as used in the present work, Box 1) or all cells by 50%, depending on the study. In recent years there has been growing recognition also of heteroresistance, describing the heterogeneity of resistance between individual clonal cells (Box 1), which appears to be very ubiquitous (4)(5)(6)(7). However, the extent to which differences in heteroresistance may contribute to MIC of a cell inoculum remains barely explored, especially relative to traditional measures such as IC 50 . The present work aimed to address this gap in our knowledge, as it could explain variability in MIC data that are used as common standards, due to parameters such as population size and the inoculum effect (Box 1).
Phenotypic heterogeneity among individual cells of genetically uniform populations is considered an evolutionarily selected strategy akin to bet hedging (8). It describes cell-tocell variation in a phenotype (e.g., stress resistance) within a cell population, resulting in phenotypically (not genotypically) distinct subpopulations among which some may be better equipped to thrive if conditions change (7)(8)(9)(10). The rationale is that generation of phenotypic heterogeneity may be a relatively low-cost strategy (metabolically) for improving survival-chances in uncertain conditions. Such heterogeneity in the case of stress resistance can be observed in almost any resistance phenotype: examples of heteroresistance are evident in microbial responses to antibiotics (2,6,11), heat (7), osmotic (10), metal (5,9), and weak organic acid (12,13) stresses. Mechanistic determinants of heteroresistance have been characterized, including cell-to-cell variation in gene expression linked to particular "high-" or "low"-variation promoters or transcription factors, commonly referred to as gene expression noise (14)(15)(16). Housekeeping genes tend to be expressed with low variation between cells, whereas high-variation promoters are more commonly associated with stress-response genes (16). Engineering increased variation of expression of an oxidative stress-response gene has been shown to increase resistance to oxidative stress in Saccharomyces cerevisiae (17), while dampening of expression variation in similar functions gave decreased resistance (9).
As food spoilage microorganisms can exhibit both high population-median resistance and heteroresistance to food preservatives (12), they provide a good example with which to address this study's main aim to characterize the relationship between MIC and heteroresistance versus population-median resistance. Among the major food preservatives, sorbic acid is a weak organic acid preservative widely used in condiments and drinks of pH #4, conditions in which the acid is undissociated and capable of traversing the plasma membrane (18). Several modes of action have been proposed for sorbic acid, including cytosolic acidification (18) and disruption of the mitochondrial electron transport chain (19). Efficacy of sorbic acid as a food preservative is hampered by the intrinsic resistance of specialized spoilage yeasts, which are capable of growth at sorbic concentrations exceeding the maximum legal limit dictated by food standards agencies (20).
Zygosaccharomyces bailii is a food spoilage yeast exhibiting a high level of resistance to weak organic acids (12). Two interspecies hybrids related to Z. bailii and also reported to spoil preserved foods are Z. parabailii and Z. pseudobailii, formed through hybridization events between ancestral Z. bailii and two unknown Zygosaccharomyces species, sharing approximately 90% to 93% sequence identity with Z. bailii (21). Z. bailii has several proposed mechanisms of resistance to weak organic acids, including a relatively thick and rigid plasma membrane (22,23), metabolism of sorbic acid to less inhibitory products (24), metabolism that is primarily fermentative (19), and a low cytosolic pH (12,20). In addition to being relatively resistant to sorbic acid (according to MIC or IC 50 ), Z. bailii also shows high sorbic acid heteroresistance (20), with a small proportion of cells capable of resisting extreme concentrations even though they may take several weeks to grow to detectable levels.
The food industry relevance of Z. bailii, coupled with its high heteroresistance makes it an ideal model for studying the relative contributions of population-median resistance and heteroresistance to MIC. In this study, we assess and compare these parameters in a panel of diverse spoilage yeasts as well as more select panels of Z. bailii, Z. parabailii, and Z. pseudobailii isolates. The work highlights the importance of context when assessing the resistance of a microorganism, with population-median resistance, heteroresistance, and inoculum size impacting the resistance of a given inoculum, nuances which are lost when an MIC at a single inoculum size is taken.

RESULTS
Use of the inoculum effect for determining heteroresistance. In order to assess the contributions made by phenotypic heterogeneity (specifically, sorbic acid heteroresistance) and IC 50 to the MIC, across large panels of species or strains, a new assay for measuring these parameters was developed based on the "inoculum effect" (25) (see Materials and Methods), to allow higher throughput than traditional measurements of heteroresistance such as from agar-based dose-response curve assays (26). Briefly, this method fitted a normal distribution curve of single-cell resistances to experimental MIC (MIC EXP ) values, by modeling MIC EXP for a given inoculum as its corresponding point on the normal distribution curve (Fig. 1B, top), e.g., MIC EXP for 100 cells/well is the point where 1 in 100 cells would survive, MIC EXP for 1,000 cells/well is the point where 1 in 1,000 cells would survive, etc. (see Materials and Methods). The median and standard deviation of the fitted normal distribution curve (Fig. 1B, bottom) could then be taken as the IC 50 and heteroresistance, respectively.
For purposes of curve fitting, lognormal distributions were adopted as these gave better fit to dose response curve data than normal distributions (Fig. S1). In addition, possible effect specifically of cell density (cell concentration) on the probability of growth occurring in a well was evaluated, e.g., due to stressor saturation or quorum sensing-type phenomena, as such effect could complicate interpretation of the inoculum effect which infers a total cell number effect. Only a minor possible effect of cell density was apparent, and similarly in high-or low-heteroresistance isolates (Fig. S2). This result supported the inference that differences measured with inoculum effect methodology primarily reflected differences in heteroresistance rather than effects attributable to altered cell density. Here, quantified as the standard deviation of a lognormal distribution of singlecell resistances. 5 IC 50 : Population-median resistance, here defined as the mean of a lognormal distribution of single-cell resistances. 6 Inoculum effect: Phenomenon by which MIC EXP increases with inoculum size, due to increased probability of resistant cells being present.

DEFINITIONS OF PARAMETERS
To further assess the use of the inoculum effect method, IC 50 ( Fig. 2A), heteroresistance ( Fig. 2B) and MIC (Fig. 2C) determinations for six strains (selected for a wide range in all three parameters based on preliminary inoculum effect measurements, data not shown) from the Zygosaccharomyces panel (Table 1) were compared between the inoculum effect method and the dose-response curve method. The parameter MIC MODEL was also taken, defined as the lowest concentration giving 99.99% inhibition of cell growth (Box 1) and consistent with previous publications in this field (12,20). This was calculated using the equation for the lognormal distribution curve fitted to MIC EXP data (as in Fig. 1B). For inoculum effect assays (carried out in broth), heteroresistance showed a moderately strong, significant correlation with heteroresistance values derived from dose-response curves (r = 0.74, P = 0.046), while IC 50 and MIC MODEL also showed strong correlations between the assays (r = 0.853, P = 0.015 and r = 0.865, P = 0.013, respectively). The inoculum effect assay was repeated on agar to assess possible effects of broth-versus agar-based assay formats (as in the above comparison) on the outcomes (Fig. S3). Similar correlations were obtained as above indicating that for all the test parameters (IC 50 , heteroresistance and MIC MODEL ), dose response results generally compare well with inoculum effect results obtained either in broth or on agar. To further examine the predictive power of this new methodology, the predicted lognormal distributions from both methodologies were compared side by side (Fig. S4). The inoculum effect distributions were relatively similar to dose response, with some outliers. As curves were infrequently directly superimposable, this methodology should be treated as a high-throughput methodology for estimating relative heteroresistance, as opposed to a precise measurement. These results supported use of the inoculum effect methodology, using a lognormal distribution for data, as an alternative heteroresistance assay to the established, but time-consuming dose-response curve method. Further assays below used the broth condition for inoculum effect assay, which offers higher throughput than on agar and is a condition also more relevant to the beverages and other product types in which sorbic acid is normally used. For stressors with a different mode of action or where there is a reasonable expectation that cell density effects are present, we would recommend testing against a dose response curve to ensure cell-density effects are not causing interference.
High-throughput determination of growth inhibition parameters for Zygosaccharomyces isolates. To investigate relationships between heteroresistance or IC 50 with MIC MODEL , we applied the inoculum effect methodology to a panel of Zygosaccharomyces sp. Isolates. To select strains for study, an initial panel of 111 yeasts isolated from foods (see Table S1) were whole-genome sequenced (WGS); most of these had previously been considered Z. bailii isolates (10,12,20). The sequencing revealed that several isolates previously characterized as Z. bailii were in fact Z. parabailii or Z. pseudobailii (see Materials and Methods). These species are known interspecies hybrids formed between Z. bailii and unknown, closely related ancestors. From these 111 isolates, a panel of 29 strains were selected for further study, to include isolates sampled from across the span of the phylogenetic tree ( Fig. S6) and diverse sources of isolation. This "Main Zygosaccharomyces panel" comprised 5, 17, and 7 isolates of Z. bailii, Z. parabailii, and Z. pseudobailii, respectively (Table 1). For parameter measurement (heteroresistance, IC 50 and MIC MODEL ) inoculum effect methodology was used, with curves fitted to MIC EXP values taken at 10 5 and 10 2 cells/well (see Table S2 for initial MIC EXP values). There was considerable variation in heteroresistance, IC 50 , and MIC MODEL among isolates of the panel, with ranges of approximately 0.03 to 0.11, 2.3 to 5.4, and 4.3 to 9.2 mM for heteroresistance, IC 50 and MIC MODEL , respectively (Table 1; Fig. 3A). Analysis for each subspecies indicated that Z. pseudobailii had a significantly greater mean heteroresistance than Z. parabailii and Z. bailii, whereas the mean IC 50 for Z. pseudobailii was significantly lower than for Z. parabailii (Fig. 3B). There were no significant differences in average MIC MODEL between the different subspecies (two sample t test assuming unequal variances), although there was both a broader spread of values and relatively low mean for Z. bailii isolates (Fig. 3B).
To further investigate the relationships between heteroresistance, IC 50 and MIC MODEL , the correlation between these parameters was calculated. Across the main Zygosaccharomyces panel of 29 strains, strains with high IC 50 tended to have high MIC MODEL and vice versa (r = 0.675, P , 0.0001), whereas heteroresistance was not significantly correlated with MIC MODEL (r = 0.326, P = 0.08) (Fig. 4A). Unexpectedly, strains with a high IC 50 tended to have a low heteroresistance (r = 0.586, P , 0.05) (Fig. 4A). When the subspecies were examined independently, Z. bailii (Fig. 4B), Z. parabailii (Fig. 4C), and Z. pseudobailii Fig. 4D) also each showed a significant correlation between IC 50 and MIC MODEL and no significant correlation between heteroresistance and MIC MODEL . In the case of Z. parabailii, it was noted that calculated correlations were strongly influenced by a single isolate (3698) with a much lower IC 50 than all the other Z. parabailii isolates (Fig. 4C). When the Z. parabailii data set was analyzed without this single isolate (n = 16), heteroresistance was strongly correlated with MIC MODEL (r = 0.660, P , 0.01) and IC 50 was not (r = 0.265, P = 0.32) (Fig. S5). The correlation between heteroresistance and MIC MODEL in this sample could suggest that heteroresistance is more strongly related to MIC MODEL when variation in IC 50 is relatively limited, as was the case in Z. parabailii after removal of isolate 3698 (see Discussion). Trends in heteroresistance, IC 50 , and MIC MODEL across a wide panel of diverse yeast species. To investigate if the trends in sorbic acid heteroresistance described above may be replicated in other spoilage yeast species, dose response curve data from a "wide panel" of spoilage yeasts isolates (Table S3) produced in a previous survey were analyzed. These isolates encompassed 26 species from 12 different yeast genera. Three Z. parabailii and two S. cerevisiae strains were present in the panel, whereas all other species were represented by one isolate only. To avoid mixing intra-with interspecies comparisons, values averaged across the three Z. parabailii and across the two  S. cerevisiae isolates were used for these two species. This wide yeast panel produced a notably wider range of IC 50 and MIC MODEL values than the main Zygosaccharomyces panel ( Fig. 5A and B). This was primarily due to comparatively low IC 50 and MIC MODEL values for certain species of the wide panel, such as Rhodotorula mucilaginosa (IC 50 ; 0.337 mM, MIC MODEL ; 0.467 mM) and Zygosaccharomyces rouxii (IC 50 ; 1.34 mM, MIC MODEL ; 2.45 mM), two species previously associated with food environments where sorbic acid is not used (27). The two panels exhibited similar variation in heteroresistance between isolates, but the Zygosaccharomyces isolates showed significantly higher mean heteroresistance (P , 0.0001) (Fig. 5B). IC 50 was very strongly correlated with MIC MODEL in the wide panel (r = 0.986, P , 0.0001) (Fig. 5C). The wide panel also showed a weaker correlation between heteroresistance and MIC MODEL (r = 0.487, P , 0.05). It should also be noted that the r value for the nonsignificant correlation between heteroresistance and IC 50 was lower than between heteroresistance and MIC MODEL , suggesting that the relationship between heteroresistance and MIC MODEL was not simply a reflection of the strong relationship between IC 50 and MIC MODEL .
The results collectively indicate a close relationship between MIC and populationmedian resistance (IC 50 ) and additionally, albeit relatively weakly, with heteroresistance, as discussed below.

DISCUSSION
This study showcases the relevance of IC 50 and heteroresistance parameters for an inhibitor's MIC against organisms, here focusing on the preservative sorbic acid. The novel, high-throughput inoculum effect methodology developed here allowed heteroresistance to be determined in a large number of microorganisms, in a manner less laborious than an established method like dose-response assay on agar. In addition, the microplate format of the new assay demands only modest physical incubator space and is amenable to automation. This methodology offers the potential for screens of heteroresistance in new ecological or industrial contexts, for example in modeling and quantitative risk assessments in the food industry. Our results indicate that a dominant indicator of differences in MIC between strains is differences in IC 50 , with heteroresistance (e.g., incidence of rare, hyper-resistant cells) playing a smaller role.
The experimental data were consistent with lognormal distributions of single-cell resistances in the yeast cell populations. Hill functions have previously been employed to model dose response curves using log 10 values for stressor concentration (10, 13), a model which assumes single-cell resistances are lognormally distributed. Here, while the fit of the lognormal distribution curve to the data was only marginally better than a normal distribution, for this study we were more interested in differences in relative distributions of single-cell resistances; therefore, lognormal was suited to the study's aims.
The possibility that cell-concentration may make a contribution to the inoculum effect that is independent of effects due to resistant-cell incidence was investigated, e.g., arising from stressor saturation or quorum sensing-type phenomena, for example. Some previous reports of the inoculum effect have not made this distinction (28,29). However, Scheler et al. (11) demonstrated that inoculum effect in E. coli for cefotaxime resulted from heteroresistance rather than cell concentration, while Steels et al. (25) ruled out certain effects of cell concentration (e.g., altered adsorption) in the inoculum effect of Z. bailii with sorbic acid. On the other hand, cell concentration affected survival probability of heat-shocked S. cerevisiae, attributed to cooperative glutathione excretion at high cell density (30). To our knowledge, no such mechanisms have been reported for sorbic acid stress. Here, analyses  Table 1). Values on the x axis indicate center of each histogram bin, where divisions between bins were at the midpoints between adjacent values shown. (B) Data adapted from (A) to compare the distributions of heteroresistance, IC 50 , and MIC MODEL among the isolates of the three subspecies. Significant differences are indicated (*, P , 0.05), using two-sample t test assuming unequal variances (n = 5, 17, and 7, for Z. bailii, Z. parabailii, and Z. pseudobailii, respectively).
indicated that cell concentration may make only a minor contribution to the increasing sorbic acid MIC MODEL as inoculum size was increased (between 10 2 and 10 4 cells per well). One possible explanation for that small effect could be related to the documented metabolism of sorbic acid by Z. bailii (24).
Heteroresistance, IC 50 , and MIC MODEL values obtained on agar (from dose response assay) showed significant correlation with the same parameters obtained in broth (from inoculum effect assay). This consistency for all three parameters whether measured in broth or agar is in line with work elsewhere reporting similar resistances to other stressors in broth versus agar (31)(32)(33). Robustness of the heteroresistance phenotype to the form of growth environment is consistent with heteroresistance being an intrinsic property of cell populations, not strongly affected by environment. This property allows resistant subpopulations to arise independent of stressor exposure, for example, so enabling prepriming for possible future perturbations. Such prepriming of heteroresistance has been described for other stressor scenarios (6,7,12), and is in line with understanding of the wider premise for phenotypic heterogeneity, or bet-hedging (8,34).
There was considerable variation in all three parameters between strains and species of the Zygosaccharomyces genus and across a wider panel of spoilage yeast species. Z. bailii is considered more heteroresistant to sorbic acid than S. cerevisiae (12). Across the main Zygosaccharomyces panel and the constituent subspecies, IC 50 was more clearly correlated with MIC MODEL than was heteroresistance. However, in the case of Z. parabailii this trend was reversed by the removal of a single outlying strain with an uncharacteristically low IC 50 value. The removal of this outlier considerably reduced the variation in IC 50 for the subspecies (standard deviation reduced to 0.055 from 0.076) but had negligible impact on variation in heteroresistance (standard deviations 0.0188 and 0.0186). This could suggest that relationships between heteroresistance and MIC are masked when IC 50 variation is high and potentially dominating effects on MIC across different organisms. Accordingly, heteroresistance may be more important for MIC in species or contexts that support lower variation in IC 50 .
Regarding the question of why one (sub-)species might evolve a higher average preservative heteroresistance than others (e.g., Z. pseudobailii, Fig. 3B), drivers could be related to the organisms' ecologies or relative incidences. This may include the size of inoculum typically associated with spoilage incidence of an organism, as larger inocula could amplify the effect that heteroresistance has on the MIC. Alternatively, the frequency of transition of the growth environment from preservative-rich to preservativepoor, as frequent changes are thought to favor high heterogeneity (34). Modeling of gene expression patterns has indicated that gene expression noise confers a stronger fitness advantage when population-averaged gene expression is further from the optimal level (35), a situation more often encountered in fluctuating environments. Regarding the particular example of this study, a less stable niche for Z. pseudobailii is partially supported by this organism being found primarily in viscous foodstuffs (Table  1), although a similar argument could be made for Z. parabailii (where average heteroresistance was lower) which colonizes solid as well as liquid foodstuffs.
There was a negative correlation between heteroresistance and IC 50 across the main Zygosaccharomyces panel. One suggestion to explain this might be that heteroresistance is more strongly selected in low-IC 50 strains in compensation for a disparity between their population-median resistance (IC 50 ) and optimal resistance when faced with stressor. That would be in line with other previously described scenarios, including the example mentioned above where noisy gene expression is especially beneficial when average gene expression is far from the optimum (35).
Other trends were evident here across a wider panel of spoilage yeast species. In particular, a weak (r = 0.487, P = 0.012) but significant positive correlation between heteroresistance and MIC emerged. This suggests that heteroresistance may make a larger contribution to MIC differences between distantly related species of spoilage yeasts than across individual isolates of the Zygosaccharomyces genus. At the same time, a correlation between IC 50 and MIC in the wide panel was stronger than in the main Zygosaccharomyces panel, possibly reflecting the wider range of IC 50 phenotypes across the wide panel. While specific trends were different between the two panels the overall conclusion was similar, that IC 50 is more closely related than heteroresistance to preservative MIC.
Why heteroresistance has not been strongly associated with sorbic acid resistance in this study requires consideration. One explanation is that if the metabolic cost to a cell of resisting sorbic acid stress is low, the fitness benefit of bet-hedging over increasing population-average resistance would be reduced. While some strategies for sorbic acid resistance such as increased expression of efflux pumps (13) would be metabolically costly, favoring heteroresistance (6), other proposed resistance mechanisms such as potential maintenance of a more fermentative metabolism (19) might confer resistance at a comparably low cost, mitigating the fitness benefit of heteroresistance over IC 50 . It is possible that stronger heteroresistance phenotypes would be found in niches where the metabolic cost of prepriming (all) cells for stress was greater. A possible related observation was overrepresentation of the high-IC 50 , low-heteroresistance species Z. parabailii in the initial Zygosaccharomyces panel used in this study (n = 71, 19 and 21 for Z. parabailii, Z. pseudobailii, and Z. bailii, respectively; Table S1). As this reflects isolation of Z. parabailii from a larger number of food spoilage events, it is consistent with the suggestion that a high IC 50 as opposed to high heteroresistance may be a more effective evolutionary strategy for growing in preserved foods, though of course we cannot discount that other traits of Z. parabailii account for its preponderance. It should also be noted that this study was limited to food spoilage isolates only and other relationships may be manifest in other ecological niches.
Conclusion. We capitalized on the inoculum effect principle to develop a new, high-throughput method for assaying heteroresistance. Application of this assay to diverse spoilage yeasts revealed that heteroresistance and especially IC 50 are related to MIC, whereas a relationship with heteroresistance was not observed when considering a panel only of Zygosaccharomyces sp. isolates. Overall, IC 50 appears to be a strong indicator of MIC, but heteroresistance may also be related to MIC between distantly related species or in contexts where IC 50 variation is low. This work highlights certain limitations of conventional MIC measurements as, due to the effect of heteroresistance, the apparent resistance of an organism could vary markedly according to inoculum size. The average resistance of a cell population, its heteroresistance, and the inoculum size should all be considered in selecting appropriate antimicrobial concentrations for specific applications, e.g., food preservation. In experimental work (e.g., in microbial challenge testing) this would translate to using an adequate population size (e.g., 10 4 cells) in the desired test volume. This parallels the relevance of heteroresistance in the clinical context, where the presence of heteroresistant-subpopulations or persister cells can lead to failure of antibiotic treatment.

MATERIALS AND METHODS
Yeast species and strains. Yeast species and strains used for experiments in this study are listed in Tables 1, S3 (alternative identifiers and origin of isolation are given where available). Members of a large, initial Zygosaccharomyces panel (Table S1) were designated Z. bailii, Z. parabailii, or Z. parabailii based on WGS analysis (see Genome sequencing) and a subset of these comprised the main Zygosaccharomyces panel used for experiments (Table 1). Several species previously designated Z. bailii were redesignated Z. parabailii or Z. pseudobailii based on the WGS analysis (below). Members of a wide panel of isolates from more diverse yeast genera (Table S3) were identified previously based on sequencing of the D1/D2 region of the 25S ribosomal DNA (rDNA) (20).
Culture conditions. The YEPD (Yeast Extract, Peptone, Dextrose) growth medium used for culturing (20 g/L d-glucose, Fisher Scientific, Germany; 20 g/L bacteriological peptone, Oxoid, UK; and 10 g/L yeast extract, Oxoid, France) was adjusted to pH 4 using 5 M HCl prior to sterilization by autoclaving. Yeasts were stored at 280°C in cryovials in YEPD medium mixed with 0.2 volumes glycerol. Cells were grown from frozen on YEPD supplemented with 2% agar (Oxoid) for 3 days at 30°C, before maintenance at 4°C for up to 3 weeks. For each biological replicate in experiments, a single colony was suspended in 1 mL YEPD broth and, after determination of OD 600 using a DS-11 FX 1 spectrophotometer/fluorometer (DeNovix Inc., DE, USA), an aliquot transferred to 10 mL YEPD broth in order to give a starting OD 600 ; 0.02. These cultures were incubated in 50 mL Erlenmeyer flasks at 24°C for .12 h with shaking (120 rev/min) prior to use in experiments at OD 600 0.1 to 2.0, when cells were growing in exponential phase. Desired concentrations of sorbic acid in the growth medium were produced using a stock solution of 20 mM potassium sorbate (Sigma-Aldrich) in YEPD, adjusted to pH 4 using 5 M HCl. Final concentrations of sorbic acid in YEPD agar or broth were produced by mixing appropriate ratios of YEPD medium (pH 4) with the 20 mM sorbic acid-supplemented medium, e.g., to make 10 mM sorbic acid in YEPD, a 50:50 ratio of YEPD medium to 20 mM sorbic acid-supplemented YEPD medium would be used. All steps were carried out under aseptic conditions.
De novo genome sequencing. De novo genome sequencing was carried out by KeyGene (Wageningen, the Netherlands) by mapping short Illumina, NextSeq2000 (San Diego, California, United States) reads to long PacBio, Sequel (Menlo Park, California, United States) reads for Z. bailii (7846), Z. parabailii (3699, containing Z. bailii-derived "A" genome and hybrid "B" genome) and Z. pseudobailii (3697, containing Z. bailii-derived "A" genome and hybrid "C" genome). Yeast isolates were grown in YEPD medium at 25°C for 48 h. Subsequently, cells were harvested by centrifugation (3,000 g, 3 min), the supernatant removed and the pelleted cultures freeze-dried prior to DNA isolation. Isolation of high molecular weight DNA suitable for long-read sequencing was carried out by KeyGene based on methods described previously (36,37).
To generate the sequencing data for Z. pseudobailii and Z. parabailii, one PacBio, Sequel SMRT cell 1M was used for each isolate. The same DNA was used for generating paired-end Illumina, NextSeq2000 reads (Flowcell P2(v2), MiSeq) that was used for polishing of the PacBio data. To obtain the initial assembly of the two hybrid species, Z. parabailii and Z. pseudobailii, the falcon-based HGAP 4 assembly method in SMRT Link portal v7.0 from Pacific Biosciences was used. As the method already includes the error correction of long reads, polishing was continued with Illumina short reads which has lower base calling error. The assembly was polished three times using Pilon software version 1.22 (38). The inverted full-size contigs in the assembly with respect to the reference genome were reverse complemented to make the assembly compatible with public reference. The contigs which were partially inverted were left as they were.
For the assembly of the Z. bailii genome, the sequencing library was prepared from PacBio high fidelity (HiFi) reads, with lower error rate (similar to Illumina short reads); hence, no polishing was required. HiFi reads were obtained via circular consensus sequencing (CSS) from the SMRT Link portal. The reads were then assembled into the draft assembly using the HiFi option of the Canu assembler. The draft assembly was polished twice with HiFi reads using Racon software (39). Then, the repetitive haplotypes were cleaned using the Purge Haplotigs pipeline (40). The remaining contigs were checked for contamination by aligning the assembly to public reference, and by aligning a part of the contigs to the BLASTn database. Finally, for all final assemblies of the three Zygosaccharomyces subspecies, a BUSCO analysis was performed to check the completeness of genomes. The pairwise alignments of genomes were done via the nucmer algorithm from the MUMmer4 package (41) and KeyGene's proprietary STL aligner (multigenome aligner). The multiple genome alignments (as pan genomes) were visualized using Mauve software (42). All sequences used in this study were submitted to the European Nucleotide Archive (ENA) under the study code PRJEB59101; de novo sequences were assigned the accession numbers GCA_949129065, GCA_949129075 and GCA_949129085 for Z. parabailii (3699), Z. bailii (7846), and Z. pseudobailii (3697), respectively.
Genome sequencing of the initial panel of Zygosaccharomyces isolates. The 111 isolates from the initial panel (Table S1) were cultured for 24 h at 25°C in YEPD and genomic DNA (gDNA) was isolated using the CTAB (cetyltrimethylammonium bromide) method (43). Some gDNA samples contained evidence of bacterial DNA, therefore, not all strains were included in the library preparation and resequencing. Resequencing was performed by Illumina HiSeq4000 with 2 Â 125 bp, on two lanes. For all Illumina sequencing, quality control was based on the QC30 threshold (80% of reads with #0.1% chance of miscalled base in all lanes of each run). The raw sequence data were trimmed on a minimum base quality of 30 (Phred scores) and sequences with one or more N-nucleotides were discarded. After trimming and filtering, read pairs for which both pairs had a minimum read length of 50 nt were kept for further analysis. The preprocessed sequencing data were mapped to the reference genome using BWA mem 0.7.17. Duplicated reads in the BAM file were marked by Picard-tools MarkDuplicates (v1.63). SAMTOOLS (v1.3) was used to view, index, sort and to merge the alignment files and the BAM files per sample genotyped using GATK 4. Variants such as SNPs and INDELs were identified using GATK4 Haplotypecaller and were stored in a single Variant Call Format (*.vcf) file. These variants were filtered on allele quality, sample quality, and allele depth. A minimum allele quality of 30, minimum sample quality of 20, and minimum allele coverage depth of 7Â were used for filtering. SNPs found in all samples that are identical were discarded as they reflect errors in the reference sequence. The filtered variants were annotated using the gene models with SNPeff. Most genes in Z. parabailii (3699) and Z. pseudobailii (3697) were assigned to the A-subgenome (highly similar to the Z. bailii 7846 genome) or to the other subgenome (i.e., B or C; see below). The long sequence aligner Blasr v1.3 was used in the Bidirectional best hit approach to identify homologs. Only coding sequences were used for this approach, extracted from the annotated gene models. Homologous genes with a high sequence similarity (i.e., above 98%) to Z. baillii genes were assigned to the A-subgenome and genes below the sequence similarity threshold were assigned to the B-subgenome (hybrid genome of Z. parabailii) or C subgenome (hybrid genome of Z. pseudobailii). Phylogenetic trees were constructed based on Genetic Distance Analysis of all open reading frames of sequences assigned to the A genome. The A genome of each sequence was compared to the A genome of Z. parabailii reference strain 3699 and an unweighted pair group method with the given arithmetic mean (UPGMA) score (44) assigned, from which a phylogenetic tree was constructed (Fig. S6). All sequences used in this study were submitted to the ENA under the study code PRJEB59101; strain-specific sample numbers are quoted in Table S1.
Dose-response curve assay for heteroresistance. Yeasts were grown to exponential phase as described above and diluted in YEPD broth to OD 600 ;0.02 (inoculum size in 75 mL ;10,000 cells) or OD 600 0.0002 (inoculum size in 75 mL ;100 cells). The use of the different inocula allowed reliable determination of % viability across a wider range of sorbic acid concentrations, with 10,000 chosen as the highest inoculum size as lower % survival values produced by plating higher inocula were not accepted by the nls() curve fitting function in R (see Fitting to normal and lognormal distributions). Aliquots (75 mL) of each inoculum size were spread across 90 mm diameter Petri dishes containing 20 mL YEPD agar supplemented with different concentrations of sorbic acid (see above). After 21 days static incubation at 23°C, colony forming units (CFU) were enumerated and percentage survival calculated by reference to control CFU counts on agar without sorbic acid. An incubation time of 21 days was chosen based on preliminary experiments of inoculum effect and dose response: 3 weeks was sufficient for outgrowth of all resistant colonies observed at week 4 (data not shown). At least three biological replicates were used in each experiment (exact numbers are clarified in figure  legends). Each of the biological replicates (from independent experiments) were the average of three technical replicates for each sorbic acid concentration. Percentage survival data were then fitted to a lognormal distribution curve in R software as detailed below (fitting to normal and lognormal distributions). Values derived with the curve fitting for IC 50 , heteroresistance and concentration at which 0.01% survival would be observed (MIC MODEL ) were recorded.
Inoculum-effect assay for heteroresistance. The inoculum effect (25) (Box 1) explains how experimentally measured MIC (MIC EXP ) can increase with size of the cell inoculum, as increased cell number increases the chance that the inoculum includes phenotypically variant cells that are hyper-resistant. Accordingly, the degree to which the inoculum effect alters MIC EXP is related to the extent of heteroresistance: MIC EXP would not be affected by inoculum size in a completely homogeneous cell population, whereas an impact of increasing inoculum size on MIC EXP will be greater the broader the range of single-cell resistances (Fig. 1A). It was reasoned that the magnitude of the inoculum effect could offer a novel, convenient measure of heteroresistance. Dose-response curves of growth-inhibitory effect typically follow a sigmoidal relationship (20) and have previously been modeled using Hill functions, (10,13), indicating a normal frequency distribution of single-cell resistances. With that assumption, a relatively small number of MIC EXP determinations at defined inoculum sizes becomes sufficient for curve-fitting (Fig. 1B, top panel), to enable determination of standard deviation (s ) from the hypothetical normal distribution of cell resistances (Fig. 1B, bottom panel). This standard deviation reflects the extent of heteroresistance, as if the difference between MIC EXP values for different inoculum sizes increases, the normal distribution will have a broader bell curve, reflecting increased heteroresistance. In addition, this curve can be extrapolated to estimate the mean of the normal distribution, reflecting the population-median resistance (IC 50 , peak of the bell curve).
Preculture of yeasts was carried out as described above. Exponential-phase cells were diluted to OD 600 ;0.2 (inoculum size in 75 mL ;100,000 cells) and serially diluted three times by 10-fold, giving a lowest OD 600 ;0.0002 (inoculum size in 75 mL ;100 cells). A total of 100,000 cells/well was chosen as largest inoculum size as this was the highest number which avoided turbidity in the well. Aliquots (75 mL) of each inoculum size were transferred to wells of a flat-bottom, 96-well microplate (CytoOne, Japan) and combined with 75 mL of YEPD broth presupplemented with sorbic acid at double the desired, final sorbic acid concentration. After sealing with insulating tape, microplates were agitated at 600 rev/min for 60 s, placed in a sealed plastic bag, and incubated statically at 23°C for 21 days (incubation time chosen as in dose response assay). The presence or absence of growth was noted visually under ;Â2 magnification using a magnifying glass. The lowest concentration at which no growth was detected in $1 of four technical replicates was scored as the experimental MIC (MIC EXP ) for that inoculum from that biological replicate of the relevant isolate.
For inoculum effect measurements with the small panel of Zygosaccharomyces strains (Table S4), the procedure was carried out as above but using inoculum sizes of 10 5 , 10 4 , 10 3 , and 10 2 cells/well, with concentration increments of 0.125 mM sorbic acid. For these particular assays, precise inoculum size was verified by plating 75 mL aliquots of the OD 600 ;0.0002 suspension (100 cells) on YEPD agar and incubating at 23°C for 7 days, the corresponding colony count was taken as the inoculum size of the lowest inoculum. The MIC EXP results were fitted to a normal or lognormal distribution curve as described below.
Fitting to normal and lognormal distributions. Royston et al. (45) describe a formula to determine the expected value of the maximum of n random variables from a normal distribution (x 1 ; ÁÁ; x n ), in terms of the theoretical mean (m) and SD (s ) of that distribution (equation 1). E½maxðx 1 ; ÁÁ; x n Þ ¼ m 1 s Ã qnormððn 2 p =8Þ=ðn 2 p=4 1 1ÞÞ (1) Assuming that the frequency distribution of antimicrobial resistance of n single cells from a homogeneous cell culture is approximated by a normal (or log-normal) one, this equation can be used to find the IC 50 (the mean) and heteroresistance (the standard deviation) using data for MIC EXP and their corresponding inoculum sizes (x). Equation 1 was fitted to experimental data using nonlinear least-squares fitting, by using the nls() function in R, to return IC 50 and heteroresistance values for each biological replicate. These were then used to calculate MIC MODEL by taking the concentration value at the 99.99th percentile of the curve (where growth inhibition would occur in 99.99% of cells; Box 1). For each parameter, the mean and standard error of the mean (SEM) across biological replicates were taken.
Plotting of curves for MIC EXP versus inoculum size was carried out using the ggplot package in R. Bootstraps to display 95% confidence interval of the line fitting were calculated using the Boot function in the 'car' package in R.
For dose-response data, sorbic acid concentration values and their corresponding percentage survival values were fitted to one minus a cumulative normal distribution function with mean and standard deviation equal to the IC 50 and heterogeneity, respectively. Fitting was done using the nls() function in R software (R x64 4.0.3). Mean and SEM of each parameter were taken from the fitted curve as for dose response curves detailed above.
Normal or lognormal distribution curves fitted to data were displayed in R (ggplot package), in the case of inoculum effect data, and were displayed using the plot() function in MATLAB R20119b in the case of dose-response curve data. Box and whisker plots, bar charts, and correlation plots were produced using GraphPad PRISM 9.0.0 software. Correlations containing points with error bars were produced using Microsoft Excel.
Simulations. Simulation of hypothetical normal distributions for single-cell resistances of organisms (for illustration of the inoculum effect concept) were produced using MATLAB R20119b. The normal probability density function (normpdf) was used to plot frequency distributions for low and high heteroresistance organisms (m = 5, s = 0.1 and m = 5, s = 0.1, respectively) with the frequency distribution divided by 10 to simulate the small inoculum size.
Data availability. All data are available in the present results section, in the supplemental material or on request. Genome sequence accession or sample numbers (ENA) are given above and in Table S1.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, DOCX file, 1.2 MB.