Introduction

The robust spatial and temporal control of symmetry-breaking and cell type differentiation is a fundamental feature of developmental systems. The mechanisms underlying the spatial control of type specialization have been extensively studied (Slack, 2014). However, there are many examples of robust cell fate choice in the absence of spatial cues. Examples include embryonic stem cells (ESCs) in culture (Pauklin and Vallier, 2013), lineage choice in the mouse blastocyst (Saiz et al., 2016; Yamanaka et al., 2010) and ‘bet hedging’ strategies in microbial populations (Maamar et al., 2007; Süel et al., 2007). In these cases, differentiation is probabilistic and is thought to depend on cell-cell heterogeneity. Intrinsic variation in cell physiology, redox state and gene expression arising from differences in cell cycle position, cell size, inheritance of cell contents upon cell division and transcriptional bursting have all been linked to fate choice (Huh and Paulsson, 2011; Raj et al., 2006). These factors can produce robust cell type proportioning if the probability of cells being in different states is predictable at the cell-population level (if the population size is sufficiently large).

Despite the widespread use of heterogeneity to facilitate developmental decision making, we still have a limited understanding of the underlying mechanisms or how variation is tuned to result in robust proportioning. Work has centred around the identification of factors that ‘prime’ cells for differentiation before they commit to a particular cell fate and express genes associated with specific cell lineages. For example, genes associated with ESC differentiation exhibit bivalent H3K4me3 and H3K27me3 epigenetic marks, suggesting a role in coordinating their expression (Azuara et al., 2006; Bernstein et al., 2006). H3K4Me3 deposition has also been shown to be sensitive to redox state, thus providing a potential explanation for the link between reactive oxygen species and ESC differentiation (Ulfig and Jakob, 2024). It has also been suggested that the level of stochasticity in the expression of lineage associated genes could influence the number of cells that are in the primed state and thus the number of differentiating cells (Desai et al., 2021). Finally, deterministic differences between cells, such as cell cycle position have been shown to influence the response to differentiation cues and lineage choice (Pauklin and Vallier, 2013).

The social amoeba Dictyostelium discoideum provides a uniquely tractable model system to study the mechanisms underlying lineage priming and the control of robust probabilistic differentiation. In response to starvation D. discoideum cells undergo a multicellular developmental programme to build a fruiting body consisting of hardy spore cells and dead stalk cells that aid dispersal from deteriorating environments. There is assumed to be selective pressure to ensure that the ‘fittest’ or most energy rich cells are chosen as spores, which requires mechanisms that bias differentiation. Furthermore, there is assumed to be strong selective pressure to ensure robust output of ‘optimal’ cell type proportioning. Excess stalk cell production is costly because it reduces the number of viable spores available for dispersal, while underproduction of stalk can compromise fruiting body architecture, limiting the success of spore dispersal (Buttery et al., 2009; Madgwick et al., 2018; Rodrigues and Gardner, 2022; Wolf et al., 2015).

Symmetry breaking and initial cell type differentiation in D. discoideum does not depend on spatial cues. Instead, cell fate choice has been associated with differences in cell cycle position at the time of starvation (Gomer and Firtel, 1987; Gruenheit et al., 2018). Cell cycle position affects the threshold of responsiveness to diffusible inducers of spore or stalk cell differentiation, such as cAMP and DIF-1 (even though all cells experience the same amount of each signal) (Chattwood et al., 2013; Gruenheit et al., 2018). Under this scenario, stochastic variation in cell cycle length generates a steady-state distribution of cell cycle positions. Optimal proportioning of stalk and spore cells is presumably reached because cell fate propensity has been evolutionarily tuned through the cell cycle when a sufficiently large number of cells is randomly sampled. It is likely that this provides an adaptive mechanism to select the least fit cells to differentiate as stalk (Zahavi et al., 2018). Cells that have just divided are primed to become stalk and are more sensitive to stalk inducing signals (and less sensitive to spore inducing signals). Post mitotic cells have just divided their cytoplasm and are likely to be energy poor compared to cells in G2. This idea is supported by the finding that glucose depletion causes cells to arrest around mitosis and differentiate as stalk cells (Gruenheit et al., 2018; Thompson and Kay, 2000).

These observations provide insights into the effects of cell cycle position heterogeneity on cell fate at a population level. However, we still have a poor understanding of how single cells interpret heterogeneity to result in probabilistic decision making or how heterogeneity is tuned to result in robust cell type proportions. To understand how cell-cell variation can be exploited to control fate choice and generate robust proportioning in the absence of cell-cell communication, we have used a combination of mathematical modelling and experimentation. The model is based on observations of probabilistic differentiation in D. discoideum. We find quantitative behaviour can only be explained by a model that not only encompasses deterministic cell cycle effects, but also the effects of stochastic cell-cell variation on the responsiveness of cells to differentiation inducing signals. Experimental observations, using single cell RNA-seq, reveals a set of genes that show extensive stochastic gene expression variation at the time of starvation. These genes show many hallmarks of lineage priming genes. They are upregulated as cells undergo a programme of differentiation and exhibit cell type specific gene expression. A genome wide screen reveals that these genes, together with cell cycle associated genes, are also more likely to be required for cell type differentiation. Stochastically expressed genes also exhibit differential H3K4 methylation. Perturbation of H3K4 methylation preferentially affects the level of their expression (and thus cell-cell variation), as well as cell type proportioning during development. These observations suggest that deterministic variation in cell cycle position acts together with stochastic variation in developmental genes to control probabilistic cell fate choice. Finally, we suggest this system is evolutionarily advantageous because it allows the use of deterministic information on the status or quality of cells (e.g., their energetic state), yet protects against potentially catastrophic effects of large-scale perturbations that cause cells to exhibit inadequate variation in deterministic properties. These studies provide further evidence supporting the key role of stochastic variation in developmental gene expression and cellular decision making.

Results

A combination of stochastic and deterministic variation explains lineage priming and fate choice in D. discoideum

To understand the mechanisms underlying lineage priming and fate choice, we first considered an earlier deterministic model that was designed to explain cell cycle dependent fate propensity in D. discoideum (Gruenheit et al., 2018). The model assumes all cells experience similar levels of a stalk-inducing factor, such as DIF-1. DIF-1 is easily diffusible and the rapid movement of cells within an aggregation result in a well-mixed population. Furthermore, levels of DIF-1 are regulated by negative feedback that can buffer perturbations in DIF-1 synthesis (Insall et al., 1992). Experimental observations instead suggest cells differ in their intrinsic threshold of responsiveness to this signal, with those above this threshold deterministically differentiating into pre-stalk cells (Chattwood and Thompson, 2011; Stevense et al., 2010). Consistent with the wealth of experimental evidence, the model assumes that these differences are partly determined by cell cycle position (Gomer and Ammann, 1996; Gomer and Firtel, 1987). It has been proposed that this is because a cell cycle associated factor (CCAF) rapidly accumulates at mitosis and increases stalk propensity (Gruenheit et al., 2018). Following mitosis, the level of CCAF decays, resulting in a decreasing propensity of cells to adopt stalk cell fate (and increasing spore propensity) as the cell cycle progresses. This model provides a relatively good fit to observed population level cell fate proportions. However, a deterministic factor governing cell fate would be expected to drive cells in the same state towards the same fate. Instead, however, cells at the same cell cycle position can adopt different fates (Gruenheit et al., 2018). Furthermore, a truly deterministic model would be expected to result in a discrete change from stalk to spore fate at the point in the cell cycle where CCAF level drops below the threshold that results in stalk fate (see Figure 1 parts A.i and A.ii). Instead, there is a gradual decrease in stalk propensity after mitosis.

Figure 1.

To develop a model that generates the sort of probabilistic differentiation observed in D. discoideum we first incorporated the deterministic influence of the cell cycle on responsiveness or sensitivity to stalk inducing factors (CCAF) (Figure 1A.i):

where A₀is the starting level of CCAF, β is its rate of decay, and t is the amount of time after the end of the previous mitosis. CCAF is measured on the scale of its biological effect rather than in terms of its molecular concentration. In theory, the decay in the level of CCAF could show a range of shapes (e.g., exponential/convex, linear, or parabolic/concave), which depends both on how the underlying factor(s) decay at the molecular level and how molecular concentrations translate into biological effects. We assume linear decay both for simplicity, and because fitting alternative models to experimental data indicate that it is the best fit model (see Supplementary text). Cells with a CCAF level above a threshold value (denoted R) adopt stalk cell fate (in response to stalk-inducing factors). As a result, stalk propensity at time t after the end of the previous mitosis (P_t) is determined by the proportion of cells with a value of CCAF (given by A_t) above R. This model would produce the steplike pattern of cell fate through the cell cycle (as described above), where all cells adopt stalk cell fate during the period of the cell cycle where CCAF levels are above the threshold (A_t ≥ R) and switch to spore cell fate when CCAF levels decay below the threshold (A_t < R)

(Figure 1A.ii). This steplike pattern will occur irrespective of how CCAF levels change through the cell cycle (i.e. regardless of whether it is linear or nonlinear) since cells in the same cell-cycle state can only be above or below the CCAF threshold that results in stalk fate.

A purely deterministic model clearly does not recapitulate the quantitative shift in cell fate propensity observed experimentally (Figure 1A.iii). We next reasoned that, for the system to show a deterministic change through the cell cycle, but result in a probabilistic output (i.e., stalk propensity), it must integrate another source of cell-cell variation that is independent of the cell cycle (along with a deterministic cell-cycle dependent factor, CCAF). We therefore hypothesised that there is also a stochastically variable cell-cycle independent factor (CCIF) contributing to sensitivity to stalk- inducing factors. The addition of such intrinsic variability makes cell fate probabilistic by allowing otherwise identical cells to adopt a distribution of possible responsiveness to stalk-inducing factors. We further hypothesized that this noisiness across cells arises from stochastic expression of a set of genes contributing to CCIF. For this, we adapt the logic of the telegraph model of gene expression (Peccoud and Ycart, 1995), where expression of gene i at time t, e_i = {0, 1}, is a Bernoulli process with probability p_i that e_i = 1 (meaning that p_i is scaled to the size of the temporal interval). Although genes might logically vary in the value of p_i, such variability does not impact our results, so we simply consider the average probability of being expressed across all genes: p̅. When gene i is expressed, it contributes x_ito the total value of CCIF, making its realised contribution g_i = e_ix_i and expected contribution in any time interval E[g_i] = p_ix_i = p̅x_i. We assume that genes vary in their contribution to CCIF according to x∼N[x̅, s²]. The total level of CCIF at time t is the sum of the g_i values for the set of G stochastically expressed genes: μ = ∑^G g_i = Gp̅x̅ and (following the law of total variance (Cornell and Benjamin, 1970)) its variance would be . This scenario results in a normal distribution of CCIF levels, CCIF∼N[μ, σ²], which we use to model the properties of the aggregate effect of all stochastically expressed genes on CCIF.

If we were to only consider the contribution of this stochastic variation, the proportion of cells adopting stalk fate (P_t) would, be determined by the proportion of the distribution that exceeds the threshold value (R) that leads to stalk fate (Figure 1B.i), which is given by the complementary cumulative distribution of CCIF values above the threshold R:

where erf is the Gauss error function, defined as . Under this scenario, the proportion of prestalk cells would remain constant throughout the cell-cycle (Figure 1B.ii), which does not follow the observed progressive shift in stalk fate probability after mitosis (Figure 1B.iii).

We integrate the deterministic influence of cell-cycle state and noisy gene expression by considering that the total level of sensitivity to stalk-inducing factors is the sum of the value of CCAF and CCIF. We achieve this integration by assuming that the mean sensitivity across cells has a constant component (CCIF) given by μ (cf. eqn. 2), and a cell-cycle dependent component (CCAF) given by A_t(cf. eqn. 1). This produces an expression for mean sensitivity at time t (C_t) that is analogous to equation (1): C_t = μ + A₀ − βt = C₀ − βt, where C₀ = μ + A₀ (i.e., is the sum of the mean CCIF value and starting CCAF value). Note that the change in average sensitivity through the cell cycle depends solely on changes in CCAF levels (βt), but the change in the average stalk propensity does not change linearly through the cell cycle because it depends on the proportion of cells with a value above R (see Figure 1.C.i). Using this expression we can rewrite the overall stalk propensity (P_t) at time t by replacing μ in equation (2) with C_t (see Figure 1C.i):

(which represents the integral over the distribution of sensitivity values above the value R at time t). This ‘stochastic-deterministic’ model generates a non-linear change in stalk propensity through time (Figure 1C.ii) that matches the probabilistic nature of cell stalk fate after mitosis (Figure 1C.iii). The model is also flexible in that, at its limits, it can produce anywhere from the step function (Figure 1A) expected for a purely deterministic process (i.e., as σ → 0) and the constant proportioning (Figure 1B) expected for a purely stochastic process (i.e., as β → 0) (see Supplementary text). To facilitate fitting of the model in equation (3) to experimental data, we can reduce the parameter space by combining the three constants (C₀, R, and σ) into a single term, denoted C^∗, which represents a sort of reference point for the model. This rescaling is achieved by assuming that the parameters are all measured in standard deviation units and that sensitivity is measured as a distance from the threshold that leads to stalk cell fate (R), such that C^∗ = C_t − R. This results in a two parameter version of equation (3), ¹ B1 + erfF(C^∗ − βt)⁄√2KJ, that can be fitted to the stalk propensity data (see below) to estimate the two parameters (C^∗and β), where the parameter estimates are in standard deviation units (meaning C^∗ represents how far mean sensitivity is from the threshold, measured as a Z-score). To test this model empirically, we started by fitting the model to existing data on stalk propensity through the cell cycle (see Supplementary Data for the original raw data and adjusted values used in the analysis). These data come from an experiment (Gruenheit et al., 2018) where time lapse microscopy was used to group cells according to their cell cycle position at the time of starvation. The fate of each cell was monitored using live cell reporter genes. For our analysis, we used the measures of relative stalk propensity over six hours (after which cells stochastically re-enter mitosis). We find that the stochastic- deterministic model provides an excellent fit to the data (adjusted R-squared = 0.92, AIC = −56.37), yielding estimates of C^∗= 0.57 and β = 0.41 (Figure 1C.iii) that correspond to an expected steady-state stalk propensity of 0.35 and an expected starting propensity of 0.72. Moreover, the stochastic deterministic model provides a significantly better fit to the data than the previously proposed model for exponential decay of CCAF (Gruenheit et al., 2018) (it is 98.35 times more likely; see Supplementary text and Supplementary Figure 1). We also evaluated the support for our model assumptions by comparing the fit of other models derived using alternative assumptions. We evaluated the assumption of linear decay in CCAF by comparison to models with exponential, quadratic, and cubic decay functions (which together capture a broad range of possible shapes), and the assumption of Gaussian variation in CCIF by fitting a model based on gamma distributed variation (which can capture a diversity of distributions, including ones approximating normality). These model comparisons all support our model assumptions (see the Supplementary text for the methods and results of these comparisons).

Cell cycle independent stochastic gene expression variation is extensive in growing cells

The stochastic-deterministic model suggests that cell fate choice in D. discoideum should depend on stochastic variation in the expression of genes associated with fate choice, as well as deterministic cell cycle dependent cell-cell variation. Single cell RNA-seq (scRNA-seq) data from D. discoideum cells prior to starvation and exposure to differentiation inducing signals has previously been used to identify cell cycle position associated gene expression variation (Gruenheit et al., 2018). Because these data are from a relatively small number of cells, but sequenced to high depth, it allowed us to determine the extent to which there is also cell cycle independent gene expression variation. We first determined the coefficient of variation (CV²) of expression for all genes. As expected, this tends to decrease as average expression level increases (Supplementary Figure 2). When this trend is accounted for, genes with greater cell-cell variability than expected for their level of expression (FDR < 0.05) can be identified (Figure 2A). Next, we used Area Under the Receiver Operating Characteristic Curve (AUROC) to identify genes with variation that can be explained because their variation in expression is associated with cells previously shown to be in early G2, late G2 or M/S phase (Supplementary Figure 2). Despite using a conservative AUROC threshold (>0.65), when these genes were removed, many genes remain that exhibit variation that cannot be explained by differences in cell cycle position (Figure 2B). Unlike cell cycle genes, principal component analysis does not result in cell groupings (Figure 2C), which suggests that this represents stochastic variation rather than a consequence of hitherto unknown extrinsic cues. Finally, we were able to determine the extent to which cell cycle associated genes exhibit stochastic gene expression. For this, the CV² of each gene was recalculated within groups of cells from each of the different cell cycle stages. Most cell cycle dependent genes were found to exhibit greater within group variation than expected (Figure 2D). This is not due to the low level of expression at specific cell cycle stages, as variation is higher at all stages, including when they are maximally expressed (Supplementary Figure 2). These results thus reveal that stochastic gene expression variation is widespread in growing cells.

Figure 2.

Stochastically expressed genes are associated with cell fate determination

To test whether genes that exhibit variability in their expression are associated with cell fate choice, we compared the developmental timing of expression of stochastically expressed genes, cell cycle associated genes and non-variable genes. The average expression of each gene during developmental stages was compared to expression during growth to generate a developmental index (where 0 is exclusive to growth and 1 is exclusive to development) (de Oliveira et al., 2019). Stochastic genes were greatly enriched (p≤0.001, binomial test) for developmental genes (index >= 0.9), whereas cell cycle and non-variable genes showed no enrichment compared to the genome wide expectation (p>0.5, binomial test) (Figure 3A). Next, we tested whether any of these groups of genes were associated with prestalk or prespore cell fate. First, RNA-seq data from prestalk and prespore cells was analysed to identify genes that exhibit cell type specific gene expression (where 0 is exclusive to prespore cells and 1 is exclusive to prestalk cells). Again, stochastic genes were significantly enriched in cell type specific genes (p≤0.001), with both prestalk (p≤0.001, binomial test) and prespore genes (p= 0.013) contributing to this enrichment (Figure 3B). Next, we determined whether these genes are more likely to be required for fate choice. We performed an unbiased large scale REMI-seq forward genetic screen (Gruenheit et al., 2021) to identify genes required for stalk cell differentiation. REMI-seq technology permits the abundance of thousands of mutants to be simultaneously quantified before or after a selection regime, such as selection imposed by the ability to undergo prestalk cell differentiation. The REMI-seq library was plated at low cell density and treated with cAMP to induce competence to differentiate followed by treatment with DIF-1 to induce prestalk cell differentiation (Figure 3C). Prestalk cells terminally differentiate as dead stalk cells after prolonged DIF-1 incubation, and thus surviving mutants with defects in prestalk cell differentiation can be enriched. After 2 and 6 rounds of growth and selection, gDNA was prepared from each biological replicate for sequencing and quantitative analysis. In order to ensure that enrichment was due to a failure to respond to DIF-1, rather than increased growth rate, we compared this set of mutants to a control selection in which cells were simply taken through an equivalent number of generations of growth (Gruenheit et al., 2021). An additional control screen was performed to identify mutants that are incompetent to differentiate at all (as either stalk or spore) because they fail to respond to cAMP. Cells were incubated in the presence of the cAMP analogue, 8-Br-cAMP, which triggers spore cell differentiation. Cells that cannot differentiate as spores were killed by detergent treatment, thus reducing their frequency in the population. REMI-seq identified 244 mutants with defects in stalk cell differentiation (Figure 3D, Supplementary Data). A subset of these mutants was recreated in the parental strain and tested in stalk cell induction assays and most mutants (8/9) exhibited defects in stalk cell differentiation (Figure 3E), illustrating the quantitative success of the REMI-seq approach. This allowed us to test whether genes required for differentiation are more likely to exhibit cell to cell gene expression variation. REMI mutants were assigned to 199 genes with intragenic insertions, or mutations where the REMI insertion lay within upstream promoter sequences (within 500bp of the transcription start site). These genes were significantly enriched for genes with variable expression (p≤0.001, binomial test) (Figure 3F). This is due to both cell cycle associated and stochastically expressed genes, as the relative number of genes identified in each class does not vary significantly from expected (chi square, p=0.77). These genes also exhibit greater variability than expected when the CV² was normalised to their expression (*** p≤0.001, t-test) (Figure 3F). These results suggest gene expression variation is a feature of genes associated with fate choice and cell type proportioning.

Figure 3.

Stochastically expressed developmental genes exhibit differential patterns of H3K4 methylation

Fate choice in D. discoideum shares features with lineage priming in embryonic stem cells, including the extensive cell-cell variation of genes associated with differentiation (Chang et al., 2008). The expression of genes associated with the differentiation is thought to be under epigenetic control due to the presence of specific patterns of epigenetic marks, including the co-occurrence of H3K4me3 and H3K27me3. The role of these modifications is not fully understood due to the difficulty with which epigenetic marks can be altered at a genome wide scale in higher organisms. The apparent absence of polycomb like proteins in D. discoideum suggests H3K27me3 modification is unlikely to play any role (Kaller et al., 2006). However, H3K4 mono, di, or tri- methylation (H3K4Me1-3), which is dependent on Set1/COMPASS, is present (Chubb et al., 2006b). Similarly, H3K9/K14 acetylation is present, which is consistent with the idea that H3K4me3 targets the Gcn5 H3K9/K14 histone acetyl-transferase to specific loci (Huang et al., 2021). Finally, recent studies have linked histone modifications to the control of transcriptional burst frequency (Weinberger et al., 2012; Wu et al., 2017), which will in turn affect the level of cell-cell variation in transcription. We therefore tested whether genes that exhibit cell-cell variation in expression in D. discoideum also show hallmarks of differential regulation by H3K4 methylation. Analysis of ChIP-seq data (Wang et al., 2021) revealed that H3K4 methylation exhibits a characteristic gene expression level dependent pattern around the transcription start site and gene body (Barski et al., 2007; Soares et al., 2017) (Figure 4A). To compare patterns in variable and non-variable genes, it was thus necessary to divide genes into bins with similar gene expression levels. Ten random samples of non-variable genes with the same distribution of expression levels as those seen in variable genes were used to compare the profile and number of genes with H3K4 methylation. This revealed variably expressed genes exhibit different profiles around the gene promoter and gene body (Figure 4B). In addition, when these were divided into stochastic and cell cycle genes, both groups were enriched for H3K4Me1 (binomial test p<0.001, Figure 4C), whilst stochastic genes were depleted for H3K4Me3 (binomial test p<0.001, Figure 4C).

Figure 4.

H3K4 methylation is required for normal expression of stochastic genes

Variably expressed genes exhibit different patterns of H3K4 methylation compared to non-variable genes. To test whether their expression also depends on H3K4 methylation, we examined a D. discoideum Set1 knockout strain in which H3K4 methylation is abolished (Chubb et al., 2006a) (Figure 5A). Bulk RNA-seq was performed on wild type and Set1 mutant cells. Most genes increased in expression, with 459 up-regulated, compared to 132 down-regulated genes in Set1^- mutant cells (Figure 5B). The affected genes were strongly enriched for variable genes identified by scRNA-seq (binomial test, p<0.01, Figure 5C). The level of enrichment increases as the stringency with which variable genes are defined is increased. This suggests that Set1 dependent genes are among the most variably expressed genes during growth (Figure 5C, Supplementary Figure 5). Furthermore, if variable genes are divided into cell cycle and non-cell cycle associated stochastic genes, Set1 dependent genes are enriched in both groups (binomial test, p<0.01) (Figure 5D), suggesting a general role in controlling the expression of variably expressed genes associated with fate choice. These findings were validated with reporter genes in which the promoters of representative Set1 dependent variable genes (GtaU, HspF-2) were used to drive RFP expression (Supplementary Figure 6). The actin15 promoter was used as a control because it exhibits little cell-cell expression variation and is unaffected by disruption of Set1 in RNA-seq. Actin15-RFP expression was similar in wild type and Set1^- mutant cells (Supplementary Figure 6), but both GtaU and HspF-2-RFP expression increased markedly in most cells (Supplementary Figure 6).

Figure 5.

Changes in gene expression can be caused by changes in burst frequency or burst size. To better understand the effects of Set1 on gene expression, we adapted the stochastic deterministic model to examine the properties of expression variation of individual genes instead of their aggregate influence on CCIF levels. For this, the probability that a gene is expressed, p_i, can be interpreted as a measure of the burst frequency, since burst frequency will dictate the probability of observing gene i being expressed in interval t. When gene i is expressed, it produces an amount of transcript given by b_i, which is a measure of burst size. Cells in the ‘off’ state are assigned a value of b_i = 0. There can be variability in the level of expression of gene i among cells in the ‘on’ state, denoted X², which can reflect both error variation in experimental estimation and biological noise/variability across cells, such that the distribution among these cells is ∼N[b_i, X²]. Hence, the average level of expression for gene i across a set of cells (which includes cells in both the on and off states) is p_ib_i, while its variance in expression would be p_iX² + p_i(1 − p_i)b_i². Consequently, an increase in burst frequency would be expected to lead to an increase in mean expression. Because both the expected mean and variance in expression levels change as a function of a change in burst frequency, we can measure the impact of a change in burst frequency on expression variability as the squared coefficient of variation (CV²), which has an expected value of (1 − p_i)/ p_i + X²/(p_ib_i²). The value of CV² would, therefore, be expected to decrease if burst frequency were to increase (i.e., Δp_i > 0) and vice versa. To test how Set1 affects gene expression, we examined the effects of Set1 disruption on gene expression in single cells. scRNA-seq data was generated from wild type and set1^- mutant cells. The level of gene expression and cell-cell variation (CV²) was determined for genes shown by bulk RNA-seq to be under the control of Set1 and by high depth scRNA-seq to exhibit significant variation in wild type cells. We find that the average expression level of these genes was affected in set1^- mutant cells (Figure 5E). As expected, based on the model predictions, the level of cell-cell variability measured by CV² significantly decreased in up-regulated genes (Figure 5F), which is consistent with an increased burst frequency. These findings are supported by quantification of GtaU and HspF-2 reporter gene expression levels and cell-cell variability (Figure 5G and H). The observed changes in gene expression are thus consistent with the idea that disruption of Set1 affects gene expression by increasing burst frequency and thus the number of expressing cells (which also reduces cell-cell gene expression variation) rather than simply further increasing the level of expression in the subpopulation of cells that already express these genes.

Loss of Set1 dependent H3K4 methylation affects cell fate choice

Genes associated with cell fate choice exhibit cell-cell variation in gene expression, they are enriched for Set1 dependent histone modifications and their normal expression is dependent on Set1. We therefore tested whether Set1 disruption affects fate choice. Disruption of Set1 did not result in defects in developmental timing, but slugs exhibit a slightly elongated and twisted morphology (Supplementary Figure 7). Clear defects were, however, evident in chimeric development between wild type and set1^- mutant cells (Figure 6A and Supplementary Figure 8). set1^- mutant cells sorted strongly towards the posterior and collar regions in chimeric slugs (Figure 6A) and to the upper and lower cup of the fruiting body, which are normally occupied by prestalk (pstB and/or pstO) cell subtypes (Supplementary Figure 8). This sorting behaviour suggests H3K4 methylation is required for normal pre-spore cell differentiation, and that mutant cells instead differentiate as pstB and/or pstO cell types. Indeed, when RNA-seq was performed on set1^- cells separated by FACS after chimeric development, set1^- cells express pre-stalk genes more strongly, whereas pre- spore genes are poorly expressed (Figure 6B). Chimeric development of wild type and set1^- strains transformed with prestalkO (ecmO-lacZ), prestalkB (ecmB-lacZ), or prespore (pspA-RFP) reporter constructs further confirm these findings. Both ecmB- lacZ and ecmO-lacZ reporter genes are more strongly expressed in set1^- cells than wild type (Figure 6C), whilst the pspA-RFP marker gene is only very weakly expressed (Figure 6C). We next used RNA-seq transcriptome profiling as it provides a sensitive method to detect defects in developmental gene expression during clonal development. Genes normally expressed in pre-stalk cells were found to be overexpressed, whereas pre-spore genes are poorly expressed (Figure 6D). Furthermore, representative pre-stalk lacZ marker genes with a long half-life, are strikingly mis-expressed throughout the slug (Figure 6E).

Figure 6.

The data from chimeric and clonal development suggest set1^- mutant cells exhibit an intrinsic defect in pre-spore cell differentiation and tendency to differentiate as pstO and pstB cell types. To determine if this is due to a lack of H3K4 methylation, rather than some hitherto unknown function of Set1, a strain was generated in which another COMPASS complex component, Ash2, was knocked out. Again, this resulted in a severe reduction in the levels of H3K4 methylation (Figure 6F). A strain was also generated in which Gcn5, which is required for H3K9/14 acetylation, was disrupted (Figure 6G). The effects of H3K4 methylation on gene expression are thought to be, in part, due to its effect on targeting H3K9/14 residues for acetylation. Indeed, disruption of set1 in D. discoideum has previously been shown to decrease H3K9/14 acetylation (Hsu et al., 2012). Both ash2^- and gcn5^- mutants phenocopied set1^- knockout cells in chimeric development with wild type cells (Figure 6G). In addition, chimeric development of a hypomorphic Set1-FLAG knock-in allele (which results in reduced di-methylation and absence of tri-methylation but little effect on mono-methylation (Supplementary Figure 9)) results in a qualitatively weaker sorting phenotype than that seen with set1^- or ash2^- mutants (Supplementary Figure 9). Together these results support the idea that H3K4 di and tri methylation is required for normal pre-spore cell differentiation.

Cell cycle position and gene expression variation interact to control cell type proportioning

Set1 disruption predominantly affects the expression of stochastically expressed developmental genes and results in defects in cell fate choice during development. This is consistent with predictions of the stochastic-deterministic model at the single cell level. A perturbation to the regulation of stochastic expression that increases the probability that genes will be expressed (i.e., increases p̅) will necessarily lead to an increase in the average level of CCIF (since the mean μ = Gp̅k̅, and such a perturbation would increase p̅), and thus will increase stalk cell differentiation. This, however, assumes that the regulation of CCAF and CCIF are independent and there are no effects on the cell cycle or the level of CCAF. To test this, we analysed the effects of perturbing cell cycle associated gene expression variation on stochastic variation, and vice versa. Cell cycle progression was blocked through cold shock, which results in widespread changes in cell cycle associated gene expression (Strasser et al., 2012). However, this did not affect the expression of representative Set1 dependent reporter genes that are not associated with the cell cycle (Figure 7A and B). Next, we perturbed the level of expression of variably expressed genes through Set1 disruption. scRNA-seq reveals set1^- cells still exhibit cell cycle associated differential gene expression with the number of M/S and G2 phase cells unaffected (Figure 7C, Supplementary Figure 10). Furthermore, cell cycle progression and timing of cell division are unaffected by Set1 disruption (Figure 7D). These data thus suggest Set1 dependent variable genes and cell cycle dependent variation are independently controlled.

Figure 7.

The stochastic deterministic model predicts an increase in CCIF should increase the likelihood of each single cell responding to stalk inducing signals (such as the stalk inducer DIF-1). To test this idea, DIF responsiveness was quantitatively compared at the single cell level in wild type and set1^- mutant strains in which a GFP reporter gene had been knocked into the DIF sensitive ecmA or ecmB prestalk genes (Supplementary Figure 11). The number of ecmB GFP positive cells was significantly higher in set1^- mutant cells at all DIF-1 concentrations tested, and the number of ecmA positive cells was higher at the lowest concentrations (Figure 8A and B). However, the level of GFP expression in each responding cell did not significantly change. This suggests that the number of cells that are in a DIF responsive state is dependent on the level of Set1 dependent stochastic gene expression (CCIF), as well as the position of in the cell cycle (CCAF). We next analysed the relationship between cell cycle phase and cell fate choice in set1^- mutant type cells. Cells were filmed growing at low density for 12-14 hours, which allowed each cell to be tracked for at least one cell division. The growth medium was then removed, and cell type differentiation induced. Cell tracking was continued, and the final fate of each cell (expression of cell type specific GFP reporter gene) was compared to the cell cycle position (timing of the last division) when differentiation was induced. Consistent with expectation from the stochastic- deterministic model (where expected stalk propensity, Gp̅x̅,increases will burst frequency), set1^- mutant cells exhibit a higher probability of adopting the pre-stalk cell fate being higher throughout the cell cycle in set1^- mutant cells (Figure 8C). Finally, we tested whether the additive effects predicted by the stochastic-deterministic model are also seen during normal development. Developmental timing and fruiting body morphology is normal when cell cycle progression is perturbed in wild type cells (by growth at 11.5°C). Similarly, although changes in the level of expression of variable genes (through mutation of set1) affect proportioning, normal fruiting bodies are formed (Figure 8D). In contrast, when both inputs are simultaneously perturbed, the effects on development are dramatic. Few fruiting bodies formed, even 24 hours after development was complete in all other samples (Figure 8D). Together, these results are consistent with a model in which symmetry breaking, fate choice and cell type proportioning in D. discoideum depends on the integration of two independent variable inputs. Since variation is widespread in biological systems, there is immense scope for natural selection to harness these properties for robust biological decision making.

Figure 8.

Discussion

The ability of isogenic cells to break symmetry and adopt different fates allows unicellular organisms to cope with dynamically changing environments and underpins the division of labour, which is a key feature of the evolution of multicellularity (West and Cooper, 2016). Cell state heterogeneity is an inevitable feature of biological systems and can provide a reliable substrate for symmetry breaking. For example, stochastic intrinsic variation can be used to predictably ‘sample’ cells from different states if population sizes are sufficiently large (even if it is impossible to know what state each cell will be in). In other cases, cells may transition reproducibly between states, such as their position in the cell cycle. Moreover, if cell state reflects properties, such as energetic or resource state, it can be used to best determine the utility of cells for different roles.

Set1 dependent methylation controls transcriptional burst frequency

Stochastic cell-cell variation can arise from noisy gene expression, which is an inevitable consequence of transcriptional bursting. Changes in burst frequency will affect the level of transcriptional cell-cell variation, as well as mean expression levels. Bursting parameters and thus levels of transcriptional noise are affected by histone modifications. Roles in differentiation are suggested by studies in embryonic stem cells, as well as other systems suggest that gene networks associated with lineage choice are often associated with specific histone modifications (Hong et al., 2011; Yadav et al., 2018). This includes the presence of marks associated with gene activation (e.g. compass complex dependent H3K4 methylation and polycomb complex dependent H3K27 methylation). Single-cell RNA-seq analysis of mouse embryonic stem cells has also shown that treatment of cells with nucleoside analogues that are removed by the base excision repair pathway result in genome wide increases in cell-to-cell variability in transcript abundance of thousands of genes (Desai et al., 2021). Importantly, this increase in transcriptional noise does not markedly change expression mean. This also affects the likelihood of differentiation, presumably because it changes the probability of cells being in a responsive state to fate inducing signals. It is, however, unknown whether this reflects endogenous control mechanisms that regulate the levels of transcriptional noise.

Our results also support the idea that gene expression noise is important for cell fate. Our studies center on Set1 dependent H3K4 methylation, which is known to recruit various HATs, HDAC’s and chromatin remodelers (e.g. ISWI, NURF, CHD) (Bian et al., 2011; Santos-Rosa et al., 2003; Shi et al., 2006; Shi et al., 2007; Sims et al., 2005; Taverna et al., 2006; Wysocka et al., 2006). Although Set1 dependent histone modifications are associated with more highly expressed loci (Barski et al., 2007; Soares et al., 2017), this relationship is not necessarily causal. In fact, genetic studies have revealed that absence of H3K4 methylation only affects a small subset of genes, and these are typically upregulated (Hsu et al., 2012; Lorenz et al., 2014; Margaritis et al., 2012; Ramakrishnan et al., 2016). In addition, previous studies have also noted a correlation between the breadth of the H3K4 methylation domain and the level of gene expression variability (Benayoun et al., 2014; Rotem et al., 2015; Sze et al., 2020; Wu et al., 2017) . Indeed, disturbance of these patterns affects transcriptional consistency (Sze et al., 2020; Wu et al., 2017). Our studies provide further support for these ideas. Moreover, they support the idea that the modification of chromatin structure provides a means by which the level of gene expression and gene expression variation can be controlled to facilitate cell fate decisions.

Noisy gene expression has the potential to buffer development against environmental perturbations

Cell fate choice and symmetry breaking in D. discoideum depend on the interplay between deterministic effects of cell cycle position and stochastically expressed developmental genes on responsiveness to fate inducing signals fate. Why should such a mechanism evolve when a purely stochastic system could achieve correct proportioning simply by random sampling of a sufficiently large population of cells? Similarly, a cell cycle system alone could also achieve correct proportioning if cells are randomly distributed through the cell cycle. D. discoideum cells that have recently completed mitosis (i.e., those that have recently divided their resources) are biased towards the stalk fate, Use of the cell cycle position presumably reflects an adaptation that helps ensure aggregations to disproportionately sacrifice ‘lower value’ cells to the stalk for building the stalk and benefit higher value (higher resource state) cells that form spores. Resource based division of labour is also seen in Sinorhizobium meliloti, where cells with highest levels poly-3-hydroxybutyrate (PHB) (which allows survival in long term starvation) do not grow and instead leave resources to be utilised by low PHB containing cells (Ratcliff and Denison, 2010).

While such examples suggest cell state heterogeneity can be utilised adaptively to make cellular decisions, it is also vulnerable to systematic perturbations that affect the proportion of cells in different states. As soil dwellers, D. discoideum cell populations are likely to be exposed to varying environmental conditions, such as temperature, pH, or available nutrients that can all potentially affect developmental signalling, cell physiology, and cell type proportions. For example, D. discoideum cells depleted in glucose or that have experienced cold shock will arrest around mitosis (and end up biased towards become stalk), which means that certain environments could lead to extreme stalk fate bias. Therefore, the evolutionary success of this organism means that mechanisms likely evolved to ensure development is buffered against environmental perturbations. Early studies identified a mechanism by which cell type proportioning could be corrected by feedback loops where the excess of one cell type leads to its trans-differentiation into the other (Belcher et al., 2022; Insall et al., 1992; Kay and Thompson, 2001). It has also been shown that disturbances that affect the cell cycle result in compensatory gene expression changes that alter the threshold of responsiveness (Chattwood et al., 2013). Our findings provide another potential solution to this problem, in which the integration of different sources of cell state variation allows systems to be robust against environmental perturbations, while still being able to adaptively exploit differences in cell state. This is because stochastic cell- cell variation, could buffer against catastrophic shifts in cell-type proportioning, by ensuring some cells always adopt different fates (Supplementary Figure 12). This would reduce the range over which feedback loops that refine proportioning need to operate. It is tempting to speculate that this complexity is a consequence of the strength of selection to ensure the correct ratio of stalk and spore cells under fluctuating natural environments that may affect signalling and cell fate choice (e.g. temp, pH, light, salt levels, humidity, food availability). Moreover, it highlights the impacts that cell-cell variation and environmental variation may have on the evolution of developmental signalling pathways. Since variation can promote or hinder developmental robustness, it is likely that many developmental systems will exhibit the goldilocks principle when it comes to adaptive exploitation of heterogeneity, where they have been evolutionarily tuned so that they are reliably and repeatedly exposed to ‘just the right amount’ of variation.

Methods

Strains, culture and development

Dictyostelium discoideum strains were derived from parental Ax3 or Ax4 strains. Cells were cultured in 10cm tissue culture dishes containing HL5 media including glucose (Formedium) or on SM agar plates (Formedium) supplemented with Klebsiella aerogenes. For development, cells in tissue culture were diluted to 1 x 10⁵ cells/ml and allowed to grow for 2-3 divisions before harvesting. Cells were induced to develop at a density of 5.1 x 10⁵ cells/cm² on non-nutrient KK2 agar (16.1mM KH2PO4, 3.7mM K2HPO4) plates containing 1.5% purified agar (Oxoid). For chimeric development, cells were mixed prior to spreading on agar. Set1 mutant lines were generated using a knock-out construct kindly provide by J. Chubb. Ash2 and Gcn5 mutants were selected from the GWDI project (https://remi-seq.org) and verified by diagnostic PCR. A strain with a random intergenic insertion (GWDI_42_D_7) was utilised as a control. For cell cycle synchronisation by cold shock, exponentially growing cells were seeded at 1.5 x 10⁶ cells/ml and incubated for 20 hours.

Reporter gene analyses

To measure cell-cell variation in gene expression of stochastic genes, promoter sequences 1kb upstream of the start codon were amplified and cloned into pDM324 (Veltman et al., 2009). The vector pDM327 was digested with NgoMIV to add sequences containing the actin15-promoter-GFP-actin8-terminator. ecmA-GFP or ecmB-KI GFP knock-in strains were generated by gene replacement (Chattwood et al., 2013). To quantify prestalk cell differentiation by FACS, growing cells were re- suspended at 1 x 10⁵ cells/ml in 1.5ml of 1x stalk solution supplemented with 5mM cAMP and incubated in 35mm plastic bottom dishes for 9hrs at 22°C. Inductions were carried out for a further 9hrs at 22°C by replenishing the cAMP and adding different doses of DIF-1. For analysis, cells were washed and collected in 1ml of KK2 + 20mM EDTA. Fluorescence was measured by FACS (Nxt Flow Cytometer (Attune)). Relative proportion of RFP and GFP positive was assessed using FlowJo software. LacZ staining of developed structures was carried out using strains expressing β- galatosidase under promoter control of fate specific markers ecmO and ecmB. Strains were fixed for 10 minutes in 4% formaldehyde and Z-buffer (10mM KCl, 60mM Na2HPO4, 40mM NaH2PO4, 1mM MgSO4). Structures were washed twice in Z-buffer before being permeabilised for 20 minutes in 0.1% NP40. Another two washes were carried out before adding staining solution (5mM K3Fe(CN)6, 5mM K4Fe(CN)6, 1mg/ml X-gal, 1mM EGTA).

Western blotting

Growing cells were lysed in nuclei buffer (40mM Tris-pH7.8, 1.5% Sucrose, 0.1mM EDTA, 6mM MgCl2, 40mM KCl, 5mM DTT, 0.4% NP40). Pellets were re-suspended at a density of 1 x 10⁷ nuclei/ml in 1x PBS including 1x protease inhibitors (Promega #G6521) and 1X SDS-PAGE buffer. Western blots were probed with antibodies to Histone-3 (ab1791 - polyclonal-rabbit IgG), H3K4-me1 (ab8895 - polyclonal-rabbit IgG), -me2 (Millipore #07-030 - polyclonal-rabbit IgG), -me3 (ab8580 - polyclonal-rabbit IgG), β-actin (Santa-Cruz Biotechnology – monoclonal- mouse IgG). Gcn5^- blots were carried out as previously described (Huang et al., 2021). Blots were imaged using a Chemidoc MP imaging system (Biorad) and Image Lab 5.2.1 software (Biorad).

Single cell analysis of cell cycle and cell fate

To follow differentiating cells, cells were transformed with the pspA-GFP cell type specific reporter gene. Log phase cells were resuspended in HL5 glucose media and 7.5 x 10³ cells were plated in 750µl deposited on a 35mm glass bottom dish (WPi). Cells were grown for a minimum of 16 hours during which time cells were imaged every 5 minutes. The time between cell divisions was determined by manually tracking single cells. After 16 hours, HL5 media was removed and cells were washed five times in KK2 buffer to remove residual growth medium. To induce differentiation, 750µl of conditioned media supplemented with 10nM DIF-1 was added. Conditioned media was collected from cells developed at a high density (1 x 10⁶ cells/ml) in 1x stalk solution (10mM MES pH6.2, 1mM CaCl2, 2mM NaCl, 10mM KCl, 200mg/ml streptomycin sulphate) supplemented with 5mM cAMP for 16 hours. Conditioned media was collected and briefly centrifuged before DIF-1 was added to washed cells. Images were taken at 5 minute intervals for another 20 hours. A final fluorescence image was taken to determine the number of prestalk and prespore cells. Cells were tracked by hand to determine the timing of the last cell division prior to the induction of cell differentiation.

REMI-seq screen for DIF-insensitive mutants

A REMI mutant library consisting of >10,000 mutants (Gruenheit et al., 2021) was grown to log phase . Cells were plated in triplicate at 2x10⁵ cells/ml in stalk medium in 10cm diameter tissue culture dishes with 5mM cAMP for 24h. Cells were then washed twice with KK2 before 24h incubation with 10nM DIF. Stalk medium was then removed and replaced with growth medium (HL5) and cells were allowed to grow until reaching confluency. Genomic DNA was prepared from the mutant library following 2 and 6 rounds of this selection and processed for sequencing. To control for mutants that are unable to differentiate in response to cAMP, an 8-Br-cAMP monolayer assay screen was performed. Cells from the mutant library were seeded in triplicate at 2x10⁵ cells/ml in stalk medium supplemented with 10mM 8-Br-cAMP. After 48 hours, detergent was added (0.1% NP40, 10mM EDTA) to remove cells that had not formed spores. Stalk medium was then removed and replaced with growth medium (HL5) and the cells were grown until reaching confluency. Genomic DNA was prepared from the mutant library following 1 round of this selection and processed for sequencing (Gruenheit et al., 2021). Analysis of mutant pools was carried out as previously reported (Gruenheit et al., 2021) using Z-score thresholds of > 1.5 for enriched mutants and < -1 for depleted mutants. Mutants (Supplementary data file 1) were identified as DIF-insensitive if they were enriched in the cAMP removal screen, not down-regulated in 8-Br-cAMP screen and not previously identified as having a growth phenotype (Gruenheit et al., 2021).

Sample preparation and analysis of bulk RNA-seq data

RNA was extracted from log phase growing cells or cells developed for 16.5 hours on 1.5% non-nutrient agar at 22°C. For RNA sequencing of chimeric development, wild type cells expressing actin15-GFP were mixed in a 1:1 ratio with set1^- cells expressing actin15-RFP and developed for 16.5 hours. Chimeric slugs were disaggregated by passing through a 25G needle before sorting (BD FACSaria) into GFP and RFP positive populations. RNA was extracted and RNA integrity number (RIN) of samples was determined by Tapestation (Agilent). Libraries were prepared from samples deemed with a RIN > 8 using an Illumina TruSeq kit and sequencing was undertaken on a Hiseq-4000 (illumina) using 100bp pair-end chemistry. Sequences were trimmed of Truseq adapters and quality controlled (Trimmomatic) by discarding reads shorter than 20bp or those where the average quality score dropped below an average of 15 in a sliding 4 base pair window. Leading and trailing bases of reads below a phred score of 30 were also removed from tags. Reads were aligned (Bowtie2) to the D. discoideum genome (an inverted repeat on chromosome 2 was masked), bamfiles were sorted (Samtools) and reads counted using the RPKM_count.py script (RSeQC). DESeq2 v1.26.0 was used for differential expression analyses. Thresholds for differential expression between samples were set at a p.adj < 0.01 with a fold change of > 2 between samples. To calculate differences in cell type specific gene expression, prestalk and prespore RNA-seq data was downloaded from the SRA (PRJNA543665) (de Oliveira et al., 2019). Genes with less than 10 read counts were removed. A cell type index was calculated for the remaining 5319 genes (Cell type index = expression count in prestalk cells / (expression count in prestalk cells + expression count in prespore cells).

ChIP-seq analysis

Bulk RNA-seq data from vegetative cells (this study) was used to rank genes based on their level of expression. Genes with detectable expression (ie > 0 normalised reads) were divided into ten equally sized bins. ChIP-seq data for two H3K4me1 and H3K4me3 replicates was downloaded from the GEO database (accession #GSE137604 Sub-Series GSE137599) as narrowPeak files ³⁹. Promoter regions around the TSS for each gene in each bin were identified (-2500bp to + 2500bp up/down-stream of the TSS) and annotated. Using functions from the chIPSeeker package and custom scripts (https://github.com/WilliamSalvidge/dictyChipSeq) annotated regions were intersected H3K4-me1 or -me3 peaks as defined by narrowPeak files. Overlaps were averaged for each expression bin and plotted using the plotAvgProf function from the chIPSeeker R package. This accounts for differing numbers of peaks in different expression bins and allows patterns of peak density between expression bins to be compared. To compare peak distribution in variable and non-variable genes, an equal number (2024) of genes was sampled using a weighted probability based on the expression of variable genes.

Identification of variable genes using single cell RNA-sequencing

Data for 81 single wild type cells isolated using the Fluidigm C1 platform were downloaded from the SRA (SAMN07833758 - SAMN0783383). Reads were normalised (DESeq2 v1.26.0) and the coefficient of variation squared (CV²) was calculated and plotted against mean expression. A trend line was fitted to the data using non-linear least squares regression (Scran v1.15.9). Genes were defined as variable (2073 genes) based on a one-sided test assuming a normal distribution around the trend but one where deviation changed depending on the mean expression of a given gene (Scran v1.15.9 - modelGeneCV2) with a FDR of < 0.05. To identify genes with a cell cycle signature, single cells were clustered using M3Drop v1.12.0 as previously reported (Gruenheit et al., 2018). Cells could be organised into three clusters of 25, 31 and 25 cells. Marker genes for each cluster were identified using logistic regression (M3Drop v1.12.0). A relaxed AUROC threshold (> 0.65) ensured that all genes possessing weak cell cycle signature could be identified (5529 genes). This allowed genes variable genes to be identified that where variation is dependent on cell cycle position (1016 genes) and independent of cell cycle position (1057 genes).

Single cell sequencing of wild type and set1^- cells

Actin15-GFP expressing wild type or set1^- cells were grown on tissue cultures dishes and harvested during log-growth. Cells were washed in 1x PBS and incubated with DAPI (2.5µg/ml). Cells were resuspended at a density of 2.8 x 10⁴ cells/ml and dispensed into a SMARTer ICELL8 3’ DE Chip using the ICELL8 cx Single Cell System. Wells of the 3’ DE Chip contain pre-printed oligonucleotides possessing well specific barcodes and UMI connected to a polydT region for hybridisation with polyadenylated transcripts. For each well of the chip, 50nl of stained cell solution was aliquoted to maximise the number of wells containing a single cell. Cells were imaged directly dispensing into nano-wells using DAPI and FITC filters and the 3’ DE Chip was frozen at -80°C. Images taken were analysed to identify wells containing individual cells based on DAPI and GFP fluorescence. In total 799 wells were identified that contained single cells (399 wild type and 400 set1^-). The 3’ DE Chip was then thawed to lyse cells and loaded onto the ICELL8 cx Single Cell System where components for reverse transcription (RT) and cDNA amplification were dispensed into chosen wells. After RT-PCR, products from separate wells were collected into a single sample, concentrated and purified according to the manufacturer’s instructions. Samples were prepared for sequencing using a Nextera XT Library Prep Kit (illumina) and sequenced on a NextSeq 500 (illumina) system utilising one flow cell and a NextSeq High-Output kit (2 x 75bp reads). One read was used to sequence the well barcode and transcript UMI, with the second reading the 3’ end of the transcript itself. This yielded >500 million reads. FASTQ files were demultiplexed (mappa), reads were then quality controlled (reads shorter than 15bp discarded, a minimum of 30% N’s allowed, phred score of 20) and trimmed of adapters (cutadapt). Reads were aligned (STAR) to the latest version of the D. discoideum genome (v2.7), sorted (Samtools) and tags counted (UMI- tools). Cells were quality controlled (Scater v1.14.6) and cells over 2 median associated deviations (MADs) from the median for library size, total number of features or mitochondrial reads excluded as outliers. This left 310 wild type and 310 set1^- cells. Genes with < 5 reads were also removed from further analyses.

Cell cycle analysis

Genes previously defined by scRNA-seq as markers of M/S and G2 cells from 81 Fluidigm-sorted cells (Gruenheit et al., 2018) were used to compare cell cycle patterns in iCell8-sorted Ax4 and Set1^- cells. The average expression of all 876 M/S genes and 642 G2 genes with expression in these samples (>0 reads) was determined. The normalised ratio of M/S to G2 marker gene expression was used to define the cell cycle position in each cell.

Model fitting to measurements of stalk cell induction

To evaluate the fit of our stochastic-deterministic model to data, we used data on stalk fate measured at different times in the progression through the cell cycle in population of cells aligned in the cell cycle (Gruenheit et al., 2018). Stalk fate propensity was measured in two genetically different sets of cells (wildtype AX3 and AX3 with a knockout of the gene gefE) grown under two conditions (‘normal’, G+, and low glucose, G−, conditions), which alter the stalk propensity of cells. We used data from the first six timepoints, corresponding to zero to five hours after the last division, and combined the four sets of cells by adjusting the values in set such that their mean propensity matched the overall global mean propensity (and hence there is no difference in average propensity of the four sets; see Supplementary Data for the raw and adjusted values). After combining the four sets, one outlier was identified (gefE⁻ under G− at hour 5), which was consistent with those cells passing a checkpoint where they re- enter mitosis, and was removed, after which data were re-normalised as described above (see Supplementary Information for a comparison of the model fitting with the outlier included). To confirm that the different sets of cells behave similarly, we also fitted the model separately to each class and see no evidence of heterogeneity of model estimates. The stochastic-deterministic model was fitted to these data using the ‘NonlinearModelFit’ function in Wolfram Mathematica version 14, which uses the Levenberg–Marquardt algorithm for least-squares curve fitting. This model fitting yielded estimates of C^∗ and β, the Akaike Information Criterion (AIC) and of the error and total sums of squares, which were used to calculate the R-squared. The details of the alternative models that were fitted and the methods used for comparing the fit of different models to the stalk propensity data are provided in the Supplementary Information.

STAR Methods

Software

• Trimmomatic - http://www.usadellab.org/cms/index.php?page=trimmomatic

• Bowtie2 - http://bowtie-bio.sourceforge.net/index.shtml

• Samtools - http://samtools.sourceforge.net

• RSeQC - http://rseqc.sourceforge.net

• Takara mappa - https://takarabiousa.github.io/mappa_userguide.html

• Cutadapt - https://github.com/marcelm/cutadapt

• STAR aligner manual - https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf

• UMI tools manual - https://github.com/CGATOxford/UMI-tools

• D discoideum genome - https://protists.ensembl.org/Dictyostelium_discoideum/Info/Index

• Mathematica 14.0 - https://www.wolfram.com/mathematica/

R - packages

• DESeq2 - v1.26.0 -

• https://bioconductor.org/packages/release/bioc/html/DESeq2.html

• Scater - v1.14.6 - https://bioconductor.org/packages/release/bioc/html/scater.html

• Scran - v1.15.9 - https://bioconductor.org/packages/release/bioc/html/scran.html

• M3Drop - v1.12.0 - https://bioconductor.org/packages/release/bioc/html/M3Drop.html

• SC3 - v1.14.0 - https://www.bioconductor.org/packages/release/bioc/html/SC3.html

• Superheat - v.0.1.0 - https://rlbarter.github.io/superheat/

Supplementary figures

Supplementary Figure 1

Supplementary Figure 2

Supplementary Figure 3

Supplementary Figure 4

Supplementary Figure 5

Supplementary Figure 6

Supplementary Figure 7

Supplementary Figure 8

Supplementary Figure 9

Supplementary Figure 10

Supplementary Figure 11

Supplementary Figure 12

Acknowledgements

This work was supported by a Wellcome Trust Investigator Award (WT095643AIA) to C.R.L.T; and grant from NERC (NE/V012002/1) to C.R.L.T and J.B.W.