Genome-level selection in tumors as a...

To disentangle the effects of therapy from the effects of natural tumor progression, we analyzed diverse datasets of both treated and untreated tumors of patients for which at least 2 samples were available, and evaluated the dN/dS ratio from the whole-exome sequencing (WES) (“Methods”). This analysis allows us to determine how treatment and natural tumor progression separately affect the selection regimes of tumor genomes, and to identify potential variations associated with treatment failure. We start with the analysis of untreated cancers, followed by analysis of treated cancers, and finally validate our findings on both published and original datasets. In the following, ‘quantifiable’ patients refer to those with at least two samples available for comparison (e.g., different primary; pre- vs. post-treatment; primary vs. metastases), and valid dN/dS estimates in both compared samples (“Methods”).

Selection in untreated cancers

To test the effect of natural progression of tumors, without therapy, we analyzed 3 published datasets of untreated cases, including 23 colorectal cancer patients⁴², 39 esophageal cancer patients⁴³ and 40 small cell lung cancer patients⁴⁵ (Fig. 1 and Figs. S1-S2). These cohorts represent a diversity of phylogenetic patterns across patients. In the colorectal cohort, early metastatic dissemination is dominant⁴², whereas in the esophageal cohort linear, punctuated, or divergent phyletic patterns are observed⁴³. All cohorts include heterogeneity among regional primary tumors.

The 21 quantifiable patients in the colorectal cohort included untreated primary and metastatic samples from each patient, and multiple primary samples for 9 patients. Metastases include both synchronous (mainly from the LN and/or liver) and asynchronous distant samples (mainly from brain and lung). We excluded post-treatment asynchronous metastases samples from 6 patients, such that the analysis represents only untreated cases. In a patient, the average dN/dS of primary samples did not differ significantly from the average value of the metastatic samples, that is, a roughly linear relationship (slope close to 1) was observed between the dN/dS values of primary tumors and metastases, even when N changed significantly (Fig. 1a and Fig. S1). A linear relationship was also observed in the comparison between early and late primary samples, that was statistically indistinguishable from the regression obtained in the progression to metastases (ANOVA, F-test) (Fig. 1a). In the esophageal cohort, multiple primary samples (at least 0.5 cm apart) were available for all 39 patients, and for 10 patients, there were also 1–3 LN samples. The analysis of this cohort also showed a wide distribution of dN/dS values across patients along with a robust linear relationship between the dN/dS values of samples in a patient, with natural progression towards LN metastases exhibiting a signature of positive selection (Fig. 1b and Fig. S2). In this case, positive selection manifests as higher dN/dS values (and typically, greater than 1) in the LN compared with the primary samples, with a steep slope, whereas the relationship between early and late primary samples is close to linear (P value < 0.018, ANOVA F-test) (Fig. 1b). In the lung cancer cohort, 3 primary tumors were collected from each patient. Comparing the dN/dS values of 36 quantifiable patients for early and late primary samples, we observed a close to linear relationship and wide distribution of dN/dS values across patients (Fig. 1c).

Overall, the analysis of these untreated cancer cases shows a wide range of dN/dS values across patients in each cancer type, with stable dN/dS values in each patient, that are largely unaffected by natural cancer progression in the primary settings. Nonetheless, natural progression toward metastases manifests cancer-specific signatures of positive (or relaxed) selection, consistent with the conclusions of previous pan-cancer analysis²⁸, and as exemplified by the esophageal cohort data, a cancer that was previously reported to exhibit signatures of positive selection⁷³.

Selection in treated resistant cancers

To identify the impact of therapy on tumor evolution, we analyzed samples before and after treatment from 5 cohorts of treated cancers, where eventually all patients developed resistance to the applied therapy. These cohorts include patients with pediatric B-acute lymphoblastic leukemia (ALL)⁴⁶, chronic lymphocytic leukemia (CLL)⁴⁸, ER-positive breast cancer³⁸, urothelial carcinoma⁵² and glioblastoma⁵⁴. The changes in N following therapy vary from cohort to cohort, yet we find a systematic effect on dN/dS as explained below (Fig. 2 and Fig. S3).

The ALL cohort included 16 (out of 20) quantifiable children who received daily standard chemotherapy (6MP; Mercaptopurine) for a year following diagnosis, and all cases recurred within several months ( <36 m). The dN/dS values before and after treatment significantly correlated, with a shift towards neutrality at relapse, whereby in cases with dN/dS < 1 before treatment, the dN/dS values tended to increase at relapse, and in cases with dN/dS > 1, there was a decreasing trend of dN/dS at relapse (Fig. 2a). This shift manifests as a statistically significantly shallower slope, compared to a linear reference model that represents the variance in the data (“Methods”) (P value < 0.057; ANOVA F-test). The post-treatment dN/dS values close to neutrality were associated with the worst outcome (rapid relapse), suggesting that the shift to neutrality is associated with higher tumor fitness and treatment failure (Fig. S4).

The CLL cohort tracks 8 patients that previously received 1 to 8 lines of standard treatment, and then were treated with venetoclax (BCL2-inhibitor) but developed resistance with median time of 15 months to relapse. Four cases developed the aggressive Richter’s transformation. Similarly to the ALL cohort, dN/dS values of 7 quantifiable patients were highly correlated with a tendency to neutral evolution following treatment at relapse (P value = 0; ANOVA F-test) (Fig. 2b). Further, patients whose pre-treatment tumors were further from neutrality (C811, C812), had the best response but eventually relapsed, concomitantly with their dN/dS values approaching neutrality. Patients with stable dN/dS values close to neutrality exhibited the worst response, whereas in cases that drifted away from neutrality (C586, C577), the time to relapse was longer (Fig. S5), consistent with the results in the ALL cohort.

The breast cancer cohort consisted of patients who developed resistance to ER-directed therapies (e.g., Aromatase inhibitors), with a signature of acquired HER2 mutations in the metastatic settings. Consistently, also in this cohort, the dN/dS values in the 7 out of 8 quantifiable patients showed a tendency to neutrality in post-treatment metastatic samples (from distant organs) compared to the primary naïve samples (P-value = 0.008; ANOVA F-test) (Fig. 2c).

The bladder cancer cohort consisted of 16 quantifiable patients who developed resistance to chemotherapy. Along with a single treatment-naive primary sample, the post-treatment samples were available from the primary tumor in some patients, whereas in others, they were collected from different distant metastases (e.g., LN, prostate, liver, lung). The dN/dS values before and after treatment were correlated, but most patients had values below 1, such that the tendency to neutrality post treatment was less pronounced than in other cohorts, although still apparent (P value = 0.098; ANOVA F-test) (Fig. 2d). Three patients received previous Bacillus Calmette Guerin (BCG) intravesical immunotherapy, which might explain some outliers.

Last, the cohort for Glioblastoma, an aggressive brain tumor with limited therapeutic options, consisted of patients who received radiotherapy and alkylating agent temozolomide (TMZ). The dN/dS values of 89 (out of 114) quantifiable patients before and after treatment were highly correlated, with a clear tendency to neutrality in the post-treatment samples (P value = 0; ANOVA F-test) (Fig. 2e). Post-treatment, many patients developed hypermutation, which was associated with a clear decrease in dN/dS and purifying selection, presumably, to mitigate the accumulation of deleterious mutations. The shift to neutrality was apparent and significant regardless of the inclusion or exclusion of these hypermutator cases. To test the association of the shift to neutrality with tumor fitness post-treatment, we performed a systematic Kaplan-Meier survival analysis, comparing patients near neutrality or approaching the neutral regime post treatment with patients that are far from or escape the neutral regime (“Methods”). The results show that the post-treatment shift to neutral evolution was significantly associated with worse prognosis, with the survival of patients improving with the distance from neutrality (Fig. 2e). These findings suggest that the neutral regime post treatment correlates with increased tumor fitness, beyond the general association of neutrality with worse prognosis in this cancer that was captured in previous analysis of primary samples from the TCGA database⁶³.

Overall, like in the untreated cases, dN/dS values in these treated cancers were largely linearly correlated. However, as therapy fails and tumors become resistant, a systematic shift towards neutral evolution was detected. This shift appears to be associated with inferior clinical outcomes, corroborating that neutral evolution corresponds to the highest tumor fitness⁵. Nonetheless, it should be stressed that the shift towards neutrality is not always apparent. For example, in cases characterized by extremely low intra-tumor heterogeneity, as previously reported for melanoma patients under chemotherapy⁵⁸, dN/dS values are expected to be highly conserved, with no apparent effect of therapy or natural progression. This is indeed the case for this cohort (Fig. S6). Thus, additional effort is required to identify the dynamics of dN/dS in a patient in such cases, as we demonstrate below.

Validation on published datasets of breast cancer patients

To validate the observations on the effects of natural progression and treatment on cancer evolution for one cancer type, we harnessed the diversity of breast cancer studies that included both treated and untreated patients, totaling 5 cohorts (“Methods” and Fig. S7). Analysis of these 5 cohorts allowed for comparisons between regional primary tumors^36,37, and between primary tumors and early local²⁹ or distant metastases^30,38, with and without therapy (Fig. 3).

The combined results from these 5 studies show that, in the untreated cases, the dN/dS values among primary tumors were strictly linearly correlated, consistent with the results on several types of untreated cancers presented above (cf. Fig. 1); however, when primary tumors were compared to metastatic samples, a shift to neutrality was revealed even in the absence of therapy (P value = 0.052; ANOVA F-test) (Fig. 3a, b). These findings were consistent with the notion of cancer-specific signatures of relaxed selection, in the natural progression to metastasis, derived from the analysis of other cohorts (see above). Considering the clinical information on the patients, the shift to neutrality appeared to originate from either early divergence of local LN in one study²⁹ or advanced stages (T4b, T4d) in another study³⁰ (Fig. 3b and Fig. S8). All the other cases, from both studies, appeared strictly linear, with only one outlier exhibiting a strong APOBEC-driven mutagenesis and high dN/dS (Fig. S8), consistent with signatures of positive selection in such cases^74,75. Further, patients under treatment who developed resistance exhibited a clear tendency to neutral evolution post-treatment, both in the primary tumors (P value = 0.044; ANOVA F-test) (Fig. 3c) and in the progression to metastatic tumors (P value = 0.008; ANOVA F-test) (Fig. 3d), in agreement with the observations on other cancers (cf. Fig. 2).

Overall, the analysis of these breast cancer cohorts presents a more complicated picture, whereby natural progression to advanced stages can also be associated with a shift toward neutrality, as a form of weakened selection, specifically, in early divergence and advanced stages, but also that treatment can cause a shift to neutrality even in the primary settings. These findings emphasize the need for further validation of the effects of therapy, for which longitudinal and temporal sampling of plasma cell cancers could be adequate¹⁵ because in these cases, it is easier to monitor the effects of therapy along the disease progression in the primary settings.

Validation on an original dataset from multiple myeloma patients

To test and validate the effect of therapy, we analyzed an original cohort of 624 multiple myeloma patients with 780 bone marrow aspirates, such that sequential biopsies across multiple stages were available from more than 100 patients (Fig. 4 and Fig. S9). The cohort spanned all stages, from very early, pre-malignant “benign” ones, including the monoclonal gammopathy of undetermined significance (MGUS) and the Smoldering multiple myeloma, followed by newly diagnosed treatment-naïve active disease (‘Pre’ state), to the early and late recurrences of the disease under continuous treatment (“Methods”). We assessed individual clinical response according to a consensus of measurements, including kappa/lambda light chains ratio (K/L) and serum M-spike from standard of care laboratory results (“Methods”). 92% of the pairs of sequential biopsies consisted of two therapy refractory (64%) or a first therapy sensitive biopsy followed by therapy resistant one (28%), highlighting the prominence of emerge of therapy resistance in MM (Fig S10). Thus, like in other treated cancers above, resistance to treatment was paramount.

Analysis of this cohort reveals dN/dS values across samples that are distributed close to neutrality, whereby disease progression correlated with an increase in the number of N mutations, MSI and CNA, and with a decrease in dN/dS (dN/dS < 1), particularly, in late relapsed refractory stages (Fig. 4a). Presumably, the decreased dN/dS is a consequence of balancing the accumulation of deleterious passenger mutations. No significant differences between stable and unstable genomes were detected in the distribution of dN/dS values (“Methods”). Notably, in the pre-treatment cases, there was a close to linear relationship between the dN/dS values, irrespective of the progression, that is, between consecutive samples of the same stage, as well as under progression from MGUS and Smoldering multiple myeloma to the treatment-naïve active disease (Fig. 4b). Comparison of pre- and post-treatment sequential biopsies confirmed a break in the correlation and a shift towards neutrality post-treatment, which persisted through early and late relapsed refractory disease (Fig. 4b). The shift towards neutrality from the natural progression in the absence of therapy to progression under therapy was apparent (P value < 0.073, ANOVA F-test).

To test the association of this shift with the clinical outcome of the patients, we performed a systematic Kaplan-Meier survival analysis with respect to the dN/dS values. Specifically, we compared all possible groups of patients in a range of dN/dS values (µ ± σ) to those outside this range. We found that deviations from neutrality post treatment were associated with better prognosis (i.e., valleys in the tumor fitness landscape), whereas under treatment the shift to neutrality was pronounced and associated with worse prognosis (i.e., hills in the tumor fitness landscape) (Fig. 4c). Additional tests further substantiated the association of this dynamics towards neutrality with poor clinical outcome. Survival tests comparing different regimes of positive against negative selection did not yield any significant correlation (Fig. S11), confirming that deviation from neutrality, either positive or negative, have equivalent positive prognostic value. Patient tumors sensitive to therapy (i.e., measurable clinical response), showed evidence of deviation from neutrality (Fig. S10), consistent with the observations on other cancers (cf. Figs. S4, S5). Cox proportional hazard analysis of the effect of different treatments on the dynamics of dN/dS confirmed that most drugs shifted dN/dS values closer to the neutral regime (Fig. 4d). Last, we performed a systematic survival analysis of paired samples, before and after treatment (“Methods”), similarly to the case of glioblastoma (cf. Fig. 2e). Patients classified as far from or escaping the neutral regime after treatment had better prognosis, whereas those near the neutral regime or approaching it after treatment had worse prognosis (Fig. 4e). This difference was more pronounced with the increasing distance from neutrality (Fig. 4e). Thus, similarly to the case of glioblastoma, disease progression under resistance to treatment correlated with neutrality and worse prognosis.

To better understand the dynamics of dN/dS throughout the evolution of tumors in individual patients, we reconstructed phylogeny, along with mutational signatures and clonal composition of each sample (“Methods”), as exemplified on several representative patients (Fig. 5a). The mutational signatures along the trees indicate that beyond COSMIC signatures 1 and 2 (that appear in most cancers), the additional significant signatures were associated with DNA damage (COSMIC signatures 3,6 and 15), and this was the case in the entire cohort (Fig. S12). Notably, dN/dS values for shared mutations at the trunk of each tree (clonal) were consistently larger than 1, clearly indicating positive selection at the onset of tumor evolution. In contrast, dN/dS values for branch mutations (subclonal enriched) were significantly below 1, indicating negative selection (presumably, of deleterious passenger mutations) in late stages. When all mutations were considered, these opposing trends largely canceled out, resulting in an effectively neutral evolutionary regime in late stages (leaves). Testing this bi-phasic tumor evolution systematically for all patients in this cohort (Fig. 5b), as well as in other cancers (Fig. S13), unequivocally demonstrated that that this pattern was universal and independent of the shape of trees or treatment status. The strong negative selection at the branches could be thought to translate into reduced clonality. However, analysis of the clonal composition indicates that the number of clones is small and highly stable (typically 2–3), with only a few cases manifesting reduced clonality (Fig. 5c). Thus, selection mostly affects the proportional size of clones, not their number. This suggests that non-genetic changes, such as epigenetic modifications, may facilitate the diversity among cells of each clone to increase (or maintain) tumor fitness, under selective pressures (“Discussion”).

Overall, these results, from individuals with sequential biopsies further validate and corroborate the expectations from theory and previous empirical evidence of the non-monotonic behavior of dN/dS during tumor progression. At initial phases (typically, characterized by low N), dN/dS increases with N due to accumulation of driver mutations, whereas at later stages (at high N), dN/dS levels off or even decreases with N due to the accumulation of deleterious passenger mutations that become a burden to the tumor and are eliminated by purifying selection^5,63,66.

Analysis of the relationship between dN/dS and tumor size

To gain further insight into the nature of the neutral regime, we analyzed the relationship between dN/dS and the respective tumor size across studies. This is because, theoretically, the neutral evolution regime following therapy could result from shrinkage of the tumors, leading to a population bottleneck and enhanced genetic drift^76,77. However, addressing this possibility by comparing the dN/dS values to the tumor size (in cases where data was available), did not provide supporting evidence (Fig. S14). Generally, there was no significant correlation between these two variables although, in several studies, neutrality was associated with large tumor size, and in only one study, we found a weak link between small tumors and neutrality (Fig. S14). Thus, the shift to neutrality appears to be associated primarily with tumor growth, rather than with population bottlenecks and drift although the latter might be a contributing factor when the size of a tumor dramatically decreases.

Relationship between dN/dS and N across patients

As described above, tumor evolution typically displays a biphasic regime in the relationship between dN/dS and N, consisting of positive correlation in initial phases, apparently, associated with the accumulation of drivers, and negative correlation in late stages, likely, due to the accumulation of deleterious passenger mutations that are eliminated by purifying selection. To further test this conclusion, we examined the relationship between N and dN/dS across the studies. In untreated cancers, dN/dS values across patients were non-monotonic in N, exhibiting an increase at low N and a decrease at high N, in primary and metastatic samples alike (Figs. S1, S2), in agreement with the expectations. By contrast, in the treated cases, this behavior was not, or at least less, apparent (Fig. S3).

The analysis of the multiple myeloma cohort provided a potential explanation for this difference. In this case, we found no clear relationship between dN/dS and N in benign and intermediate stages of the disease, but a non-monotonic behavior in N was evident in disease states (Fig. S9). Furthermore, in samples under continuous treatment, this non-monotonicity was perturbed, presumably, because treatment drives tumor evolution toward the neutral regime (Fig. S9). Thus, the present analysis corroborates on single cohorts the expected non-monotonic behavior of dN/dS in N and additionally suggests that therapy can perturb this behavior as it attracts tumors to the neutral regime. We next examined and validated these conclusions on individual patients.

Analysis of the relationship between dN/dS and N in individual patients

To further validate the dynamics of dN/dS and demonstrate how it can inform our understanding of the progression of tumors in individual patients and their response to treatment, we analyzed several case studies of patients with multiple samples. We focused on patients with at least 5 samples, where we could evaluate the progression from the treatment-naive state to the metastatic state post-treatment, along with the correlations between dN/dS and N as function of allele frequency (AF) of mutations. Systematic exploration of different AF ranges provides an approximate time signature that increasingly shifts from early, clonal-enriched mutations to late, subclonal-enriched mutations (Fig. 6a). The analysis of these patients provides further evidence of the general non-monotonic character of tumor evolution, whereby signatures of positive correlation between dN/dS and N were mostly observed in early stages, whereas signatures of negative correlation were observed in late stages; in addition, however, patient-specific behaviors were demonstrated (Fig. 6a).

Specifically, we found support of the non-monotonic tumor evolution regime in many patients. For example, such evidence was obtained for two patients with the largest number of samples from the colorectal cancer cohort⁴² that included both treated and untreated samples: patient V930, representing the common case of early metastatic dissemination, and patient V750, representing the rarer case of late dissemination, as deduced by the original phylogenetic analysis⁴² (Fig. 6b and Fig. S15). In both cases, the dN/dS values increased toward neutrality as the tumor progressed, with the treated metastatic samples being the closest to neutral, corroborating the attraction to the neutral regime in an individual patient. The attraction to the neutral regime as the tumor progressed under treatment blurred to some extent the correlations between dN/dS and N, for the treated samples; however, for the untreated samples, both the positive and the negative correlations, at early and late stages, respectively, were evident and stable (Fig. 6b and Fig. S15).

Although these observations appeared to be robust with respect to the evolutionary path (early vs. late dissemination), the impact of treatment and the timing of driver mutations acquisition varied between cases. For example, in the patient with the largest number of samples from a bladder cancer cohort⁵², the effect of chemotherapy on the clonal dynamics was substantial (Fig. 6a). In this patient (P117), there was a significant reduction in N following treatment, in both primary and metastatic lesions, compared to the single treatment-naïve primary sample. However, dN/dS in both the untreated and treated primary samples was close to neutrality, with somewhat lower values in the treated samples (Fig. 6c). After treatment, N increased significantly from the primary samples to the metastatic ones, with a signature of purifying selection that was consistent with the high mutational burden and manifested as lower dN/dS values in the metastatic samples (Fig. 6c). The drastic effect of therapy makes it harder to identify the accumulation of driver mutations and the post-treatment attraction to neutral regime. Nonetheless, the systematic examination of the relationship between dN/dS and N identified a positive correlation between dN/dS and N that was dominant in relatively late stages (Fig. 6a, c). Thus, despite the overall reduced dN/dS values, this dominant positive correlation in these later stages suggests that beyond the advent of the early drivers, accumulation and selection of late (and likely resistant subclonal) drivers continue to occur and shape the evolution of the tumor genome under treatment, consistent with the conclusions of the original study⁵².

Last, as another demonstration of how the dN/dS analysis can help assessing the response to treatment, we analyzed a melanoma patient that was tracked over 9 years, 6 years under immune blockage inhibition, including 37 samples, from the treatment-naïve primary stages, through treatment, to the last metastatic burst⁵⁹ (Fig. 6a). The original phylogenetic analysis of the samples identified 6 lineages, from one of which the last metastatic burst originated. Analysis of the data from this patient showed that, despite the general low dN/dS values (expected for melanoma due to the typical high mutational burden) and the general low intra-tumor heterogeneity⁵⁸ (cf. Fig. S6), the lineage that resisted therapy originated from those samples that most closely approached the neutral regime (Fig. S16). Thus, this case study provides additional evidence of how relaxed selection is indicative of tumor progression under treatment, even in cases with a narrow dynamic range of dN/dS; furthermore, phases of positive selection in tumor evolution were identified (Fig. 6a and Fig. S16).