Measures of genetic diversification in somatic tissues at bulk and single-cell resolution

Abstract

Intra-tissue genetic heterogeneity is universal to both healthy and cancerous tissues. It emerges from the stochastic accumulation of somatic mutations throughout development and homeostasis. By combining population genetics theory and genomic information, genetic heterogeneity can be exploited to infer tissue organization and dynamics in vivo. However, many basic quantities, for example the dynamics of tissue-specific stem cells remain difficult to quantify precisely. Here, we show that single-cell and bulk sequencing data inform on different aspects of the underlying stochastic processes. Bulk-derived variant allele frequency spectra (VAF) show transitions from growing to constant stem cell populations with age in samples of healthy esophagus epithelium. Single-cell mutational burden distributions allow a sample size independent measure of mutation and proliferation rates. Mutation rates in adult hematopietic stem cells are higher compared to inferences during development, suggesting additional proliferation-independent effects. Furthermore, single-cell derived VAF spectra contain information on the number of tissue-specific stem cells. In hematopiesis, we find approximately 2 × 10⁵ HSCs, if all stem cells divide symmetrically. However, the single-cell mutational burden distribution is over-dispersed compared to a model of Poisson distributed random mutations. A time-associated model of mutation accumulation with a constant rate alone cannot generate such a pattern. At least one additional source of stochasticity would be needed. Possible candidates for these processes may be occasional bursts of stem cell divisions, potentially in response to injury, or non-constant mutation rates either through environmental exposures or cell-intrinsic variation.

Keywords: evolutionary biology; evolutionary inferences; healthy human tissues; human; sampling; single-cell mutation burden; stem cell dynamics; varient allele frequency.

Figure 1.. The distribution of variant allele frequencies changes with the growth phases and by sampling.

(a) In the current population, cells divide symmetrically into two daughter cells or asymmetrically with only one daughter cell kept in the focused population. All other events are mathematically equivalent to and are treated as a part of cell death. (b) The rates of symmetric and asymmetric division change during the population growth and lead to a dynamic distribution of variant allele frequencies. (c) The observed VAF distribution is shifted again during sampling compared to the VAF of the whole population – a fact should be considered when inferring population properties through genetic data.

Figure 2.. Bulk sequencing based variant allele frequency (VAF) and mutation rate inferences in healthy esophagus.

(a) Expected VAF distributions from evolving Equation 1 to different time points for a population with an initial exponential growth phase and subsequent constant population phase (mature size $N\mathrm{=}{\mathrm{10}}^{\mathrm{3}}$). Once the population reaches the maximum carrying capacity, the distribution moves from a $\mathrm{1}\mathrm{/}{f}^{\mathrm{2}}$ growing population shape (purple) to a $\mathrm{1}\mathrm{/}f$ constant population shape (green). Note that the shift slows considerably at older age. (b) VAF from healthy tissue in the esophagus of nine individuals sorted into age brackets. The youngest bracket, 20–39, is closer to the developmental $\mathrm{1}\mathrm{/}{f}^{\mathrm{2}}$ scaling. The older age brackets are both close to the constant population $\mathrm{1}\mathrm{/}f$ scaling, resembling the theoretical expectations. (c) Inferred mutation rates increase linearly with age. (d) Simulations of slowly growing stem cell populations reveal that mutation rates appear to increase with age, although the true underlying per division mutation rate remaining constant (see Figure 2—figure supplement 1 as well).

Figure 2—figure supplement 1.. Comparison of an exponential and logistic growth model.

(a) Population size over time for exponential (dotted) and logistic (dashed) growth functions. At the point of maturity ($t\mathrm{=}\mathrm{20}$) the exponential function reaches the population capacity ($N_M\mathrm{=}\mathrm{10}\mathrm{}\mathrm{’}\mathrm{}\mathrm{000}$), and the logistic function equals $\mathrm{0.99}\mathrm{}{N}_M$. The time points at which the VAF spectra are measured are indicated by solid vertical lines: $t_{\mathrm{1}}\mathrm{=}\mathrm{25}$ (blue) and $t_{\mathrm{2}}\mathrm{=}\mathrm{50}$ (orange). (b) Comparison of the expected variant allele frequency (VAF) spectra, calculated with Equation 1, in the exponential (dotted) and logistic (dashed) growth models measured at the time points $t_{\mathrm{1}}$ and $t_{\mathrm{2}}$. For reference, the theoretical predictions for an exponentially growing population without cell death and a constant population in the long time limit are shown as solid purple and green lines, respectively.

Figure 3.. Inference of evolutionary parameters on simulated stem cell populations.

Simulated populations were run up to age 59, growing exponentially from a single-cell to constant size $N_M={10}^{\prime }000$ at age $t_M=5$, with mutation rate $\mu =1.2$ and division rates $\lambda =5$ and $p=0.4$. Where sampling is mentioned, the sample size 89 was taken. (a) The single-cell mutational burden distribution. The compound Poisson distribution (dashed line) matches the burden distribution when averaging over multiple independently evolved populations (filled curve). (b) Distribution of estimated mutation rates from 10’000 individual simulations, obtained from burden distributions of the complete populations (blue) as well as sampled sets of cells (orange). Because the expected mutational burden distribution is unaltered by sampling, the expected estimate of the mutation rate from Equation 5 remains unchanged: $E({\tilde{\mu }}_{pop})=E({\tilde{\mu }}_{sample})$. However, sampling increases the noise on the observed burden distribution, which results in a higher error margin of the estimate: ${\displaystyle \sigma ({\tilde{\mu }}_{pop})<\sigma ({\tilde{\mu }}_{sample})}$. (c) VAF spectra measured in the complete population (blue) and a sampled set of cells (orange). In contrast with the mutational burden distribution, strong sampling changes the shape of the expected distribution. A single simulation result is shown (diamonds) alongside the theoretically predicted expected values for both the total and sampled populations (Equations 1 and 12) (dashed line) and the average across 100 simulations (solid line). (d) Distribution of $N_M$ and $p$ inference results for 100 simulated and sampled populations, through estimation of $\tilde{\mu }$ and $\tilde{\lambda }$ from the single-cell burden distribution and fitting the number of lowest frequency ($1/S$) mutations to the theoretical prediction in Equation 1 (see Figure 3—figure supplements 1–3 as well).

Figure 3—figure supplement 1.. Stochastic simulations of the Variance/Mean of the mutational burden distribution over time for a per cell division mutation rate of μ=1.3 and varying stem cell population size N and asymmetric division probability p.

Stochastic fluctuations are pronounced for small population size $N$ and low asymmetric division probability $p$.

Figure 3—figure supplement 2.. Likelihood of the Variance/Mean to be in the interval 3<θ<5 for a per division mutation rate of μ=1.3.

Figure 3—figure supplement 3.. Likelihood of the Variance/Mean to be in the interval 3<θ<5 for a per division mutation rate of μ=3.

Figure 4.. Evolutionary inferences in single-cell hematopoietic stem cell (HSC) data.

(a) The single-cell mutational burden distribution of the data (bars) and the compound Poisson distribution obtained from its mean and variance, used to obtain the estimated per division mutation rate $\stackrel{\mathrm{~}}{\mu }$. (b) Distribution of mutation frequencies of the data and theoretically predicted average fitted to only the lowest frequency ($\mathrm{1}\mathrm{/}S$) data point. (c) Difference $\mathrm{\Delta }\mathrm{}{v}_f$ between the measured value of the VAF spectrum at the lowest frequency ($\mathrm{1}\mathrm{/}S$) and its prediction from Equation 1, for varying total population size $N$ and asymmetric division proportion $p$, with fixed maturation time $t_M\mathrm{=}\mathrm{5}$ and operational hematopoietic population size ${\displaystyle N_H=50}$. The solid line denotes the plane of best fit where this difference is 0. (d) Maximally inferred population size $N$ (taking $p\mathrm{=}\mathrm{0}$ in (c)) for variation of the maturation time $t_N$ and the operational hematopoietic population size $N_H$ (see Figure 4—figure supplement 1 as well).

Figure 4—figure supplement 1.. The standard deviation on the variant allele frequency (VAF) spectrum increases for higher frequencies.

(a) The VAF spectrum $V_f$ averaged across 100 simulations of a population evolved according to the model described in section 1.1. The standard deviation from the mean is shown as the band $V_f\mathrm{\pm }{\sigma }_f$ around the average. (b) The standard deviation across all simulations for each frequency $f\mathrm{=}\mathrm{1}\mathrm{/}N\mathrm{,}\mathrm{2}\mathrm{/}N\mathrm{,}\mathrm{\dots }$, scaled by the average spectrum $V_f$.