Interplay between whole genome doubling and the accumulation of deleterious alterations in cancer evolution

LOH and WGD are common events in NSCLC

First, we explored WGD, LOH and mutational burden in NSCLCs from the TRACERx¹⁶ and TCGA datasets¹⁷. WGD was a frequent event in NSCLC (64.45% in LUAD, 66.12% in LUSC tumors), consistent with previous work^1,16. Using TRACERx data, where multiregion data were available, we confirmed that WGD is primarily an early event in cancer progression, with only a small subclonal proportion (6.25% in LUSC and 1.63% in LUAD) as previously reported in Jamal-Hanjani et al. (2017) (Figure 1a,b).

We noted a lower proportion of the genome subject to LOH in non-WGD tumors compared to their genome doubled counterparts (Figure 1a-c), but as expected, by virtue of the additional genome copies, on average only 1.32% of the genome was haploid in WGD tumors, compared to 9.55% in nWGD cases (Figure 1a,b). Clonal LOH was significantly more common than subclonal LOH (average of 25.85% vs 5.96%, Figure 1c), suggesting the majority of LOH occurs early, likely before WGD.

Additionally, we observed a significant positive correlation between the amount of haploid LOH in non-WGD tumors and the frequency of WGD events across tumor types (Figure 1d). These observations led us to explore whether WGD could mitigate the potential negative effects of haploidy and the accumulation of deleterious alterations (mutations and/or SCNA) in these regions.

The fitness cost of WGD is offset by its positive impact

To explore the impact of a WGD event on cancer cell viability and fitness and whether a duplication event may be selected to buffer deleterious alterations (Figure 2a,b), we adapted a tumor progression model from¹⁴.

In our model, alterations in cancer genes is associated with a fitness gain (s_d) while passenger alterations in haploid regions are associated with a weak fitness cost (s_p), reducing the birth rate. WGD influences the fitness cost of these passenger alterations, reducing their negative impact, but is itself associated with a fitness cost (s_WGD).

To explore WGD we simulated the evolution of cancer populations (starting with a 1,000 cells), and varied the fitness costs of passenger alterations (s_p ∈{0-0.01}, where 0.01 represents a 1% increase in waiting time to birth) and WGD (s_WGD ∈{0-1.5}) (Supplementary Table 1). Deleterious passenger alterations represent either mutations or SCNA losses affecting genes in regions of LOH. To explore whether WGD is selected, at the end of each simulation we calculated the proportion of WGD in the population. We initially assume the rate of deleterious alterations, u_p, is 0.5 and that cells adapt to the cost of WGD after ten generations. The deleterious alteration rate reflects the product of number of putatively deleterious sites (i.e. the extent of haploid LOH before doubling) and the alteration rate per-base-per-cell-division. We simulate until the population reaches 20,000 cells or 2,000 generations. Using these parameters, after 1,000 generations we observed an average of 489 [381-514] passenger and 8 [7-9] driver alterations.

Our simulations suggest a relationship between fitness cost of passenger alterations and the likelihood of WGD being selected; higher fitness cost of passenger alterations was associated with increased WGD (Figure 2c). However, this relationship was also dependent on the cost of WGD; a higher fitness cost was associated with reduced WGD, when passenger alterations were associated with very minor fitness costs (s_p <0.00002). Notably, provided the fitness cost of passenger alterations exceeded 0.0002, WGD was selected, even when the fitness cost of WGD was high (s_WGD>0.5). Moreover, if the fitness cost of WGD was low (s_WGD=0.05) and the cost of passengers relatively high (s_p=0.001), WGD was selected even when the fitness cost was maintained for >100 generations (Supplementary Figure 1a).

However, when the average passenger fitness cost was negligible (sp <0.00001), the benefit of WGD did not compensate for its cost. Likewise, if the deleterious alteration acquisition rate was low (u_p<0.2 for s_p = 9.1 × 10^-5 or u_p<0.05 for s_p= 4 × 10^-4), e.g. due to a low proportion of the genome subject to LOH, or if the WGD fitness cost was high (sWGD>0.5) and durable (>25 generations) WGD was rarely selected (Figure 2d).

Notably, a stable genome without any haploid regions, will exhibit a deleterious acquisition rate approaching zero (i.e u_p ≃ 0). This may be the case for hyper-mutator tumors (for example POLE mutant tumors), which display a lower burden of SCNA and low frequency of WGD¹⁸. By contrast, a high deleterious rate of deleterious alteration acquisition may occur in with a high degree of haploidy and chromosomal instability, whereby there is a high probability of a SCNA leading to homozygous loss. These results are in agreement with the relationship between LOH and WGD proportions across cancers (Figure 1d).

To evaluate whether the fitness costs explored are compatible with previously documented lack of detectable negative selection¹⁹ we explored the relationship between the cost of passenger alterations and purifying selection. Despite being associated with WGD, weakly deleterious alterations (s_p<0.0005) largely escaped detectable negative selection (Supplementary Figure 1b). Conversely, alterations with a substantial fitness cost (s_p >0.01) were associated with strong negative selection.

Taken together, this model, despite simplifying cancer evolution, suggest WGD in cancer cells can act as an evolutionary mechanism that might serve to buffer the deleterious effect of somatic alterations. The model suggests that provided there is a sufficient deleterious alteration rate, through SCNA or mutations, WGD can potentially be selected.

Accurate timing of mutations pre and post WGD

To quantify the impact of WGD in cancer evolution using sequencing data and to determine whether there is a shift in detectable selection following a duplication event, it is imperative to be able to time somatic alterations in relation to WGD. Previous work has suggested WGD provides a natural mechanism to temporally dissect mutations^20,21. In brief, mutations occurring prior to WGD should be present at multiple copies, while those occurring after a doubling event would only be present at one copy (Figure 3a).

However, this assumption has not been subject to experimental validation. Therefore, we utilized an isogenic genome doubling model system involving genome-doubled HTC-116 clones deriving from a non-genome doubled common ancestor³ and exome-sequenced the ancestor in addition to two diploid and four tetraploid cell lines at two time points (passages 4 and 50). Given that the common ancestor was diploid, all clonal mutations should occur before WGD, while private mutations should almost all occur after WGD.

Reassuringly, 90.19% of pre-genome doubled mutations were correctly timed as such. Conversely, 94.70% tetraploid private mutations were correctly classified as occurring after WGD (Figure 3b). Applying temporal dissection of mutations to tumors exhibiting WGD in TRACERx and TCGA datasets suggests that the majority of detectable clonal mutations in NSCLC accumulate prior to WGD (Figure 3c).

Purifying selection on essential genes prior to duplication

If WGD has a significant impact upon selection and the evolutionary course of the disease, we reasoned one would expect to see a difference in the selection of mutations before and after doubling. While our simulations suggest weakly deleterious alterations may escape detectable selection, purifying selection will likely operate on clones bearing deleterious alterations in housekeeping or essential genes, located in haploid regions of the genome (where s_p > 0.01). We reasoned that if WGD buffers deleterious alterations, purifying selection pressures will be relieved after duplication of essential genes.

To investigate the selective pressures acting before and after WGD, and the effect of LOH, we applied a modified dNdS ratio test¹⁹ to early (pre-WGD) and late (post-WGD) mutations within segments of LOH. Under the assumption that synonymous mutations are neutral, this ratio can be informative about the direction of selection: ratios >1 indicate positive selection, while ratios <1 are consistent with negative selection. We investigated different sets of genes, including essential genes identified by using mutagenesis screens in haploid human cells²² and LUSC and LUAD-specific cancer genes described in the literature^19,23–26. We focussed on dNdS values for truncating mutations (including nonsense and splice-site mutations), which would be most closely associated with protein dysfunction and should in theory be subject to stronger purifying selection (Figure 4a). dNdS values for missense mutations are shown in Supplementary Figure 2.

Analyses at the whole exome level (all genes) showed no deviations from neutrality (Figure 4a), as previously reported¹⁹, and consistent with our simulations (Supplementary Figure 1). However, signals of purifying selection were observed for essential genes occurring in early mutations (pre-WGD) encoded within genomic segments of LOH (Figure 4a, Supplementary Figure 3). These data suggest homozygous disruption to essential genes in regions of LOH in diploid genomes resulting in disruption of both alleles, leads to a detectable fitness impairment; indeed, on average, approximately half of clones harboring early truncating mutations in regions of LOH are predicted to be lost due to purifying selection. However, due to the limited number of point mutations in essential genes, we can estimate that on average only one clone harboring a nonsense mutation in an essential gene is eliminated per tumor prior to WGD during its evolution.

We hypothesized that following a WGD event, new mutations in genomic segments of LOH occur less frequently in a haploid context and therefore purifying selection would be weaker. Consistent with this, mutations occurring after WGD in essential genes were not found to be subject to significant purifying selection (Figure 4a). Likewise, in non-LOH regions (which harbor more mutations), dNdS values were not significantly different from 1. These results suggest that genomic segments that accumulate truncating or nonsense mutations but still maintain one or more wild-type copies of the gene are not subject to strong negative selection as they are still viable. Therefore, by reducing the haploid genome fraction of cancer cells, WGD may act to mitigate the detrimental effect of deleterious mutations.

To rule out an excess of SNP contamination in essential genes leading to spurious signatures of negative selection, we compared the proportion of somatic mutations contained in the SNP database (dbSNP) in essential²² vs non-essential genes and all genes (Supplementary Figure 4). The proportion of mutations classified as SNPs was low in essential, non-essential and all genes in both LUAD (0.017, 0.022 and 0.022, respectively) and LUSC (0.025, 0.026 and 0.025, respectively), suggesting germline contamination is not responsible for the observed negative selection.

In order to assess the effect of using different essential gene lists, in addition to the essential genes from²²(Figure 4a), we took advantage of CRISPR screens and calculated the dNdS ratio for pre-WGD mutations in LOH in NSCLC specific essential genes²⁷. We noted strong purifying selection acting on genes with the highest gene-knockout effects (i.e. genes whose deletion has the most negative impact on cell viability) (Supplementary Figure 5). We also quantified selection in other essential gene lists from additional CRISPR-based model-systems^28,29. The results were broadly consistent. For TCGA-LUSC, we observed the same trend regardless of the gene-set used (Supplementary Figure 6). For TCGA-LUAD we noted that while negative selection was observed prior to doubling in all cases, for genes derived from²⁹, a weak signal of negative selection was observed for post WGD mutations, potentially due to some genes acting in a haploinsufficient manner.

Taken together these results suggest that WGD may mitigate against nonsense mutations in essential genes. However, given that these mutations are relatively rare events during cancer evolution (u_p<0.0001), these data, together with the simulations, suggest WGD cannot be selected solely to buffer against nonsense mutations in essential genes.

We therefore next explored whether WGD may also have an impact on the acquisition of SCNAs. Homozygous deletions are rare in NSCLC genomes, consistent with complete loss of gene function generally having a negative effect on cell fitness. Genomic segments which are haploid may be particularly susceptible to copy number loss events, which would result in homozygous deletions. To explore whether the SCNA landscape is influenced by negative selection we considered whether the proportion of haploid LOH in WGD tumors was less than expected by chance. We simulated copy number losses randomly (at the chromosome arm level) in WGD tumors, fixing the number of chromosome arms in pre-WGD LOH based on the observed data. Next, we quantified the number of post-WGD losses resulting in haploid LOH. Through simulations in WGD tumors exploring where copy number losses occur across the genome, we found significantly less haploid LOH than would be expected by chance, consistent with negative selection removing those tumor clones (Figure 4b). Furthermore, we explored the correlation between copy number losses and essentiality of the chromosome arms (as captured by the density of essential genes)³⁰ and observed a significant negative correlation between the essential score and haploid LOH occurrence for each chromosome arm (Figure 4c), consistent with negative selection resulting in fewer copy number losses where there is an abundance of essential genes. This correlation, however, was not observed when we focused on copy number losses occurring after the WGD event (Figure 4d), consistent with negative selection sculpting the SCNA landscape, reducing haploidy and homozygous deletions.

Similar selection patterns across other cancer types

Prompted by these findings, we sought to further investigate selection across other cancer types. Consistent patterns were observed among most of those cancer types with sufficient patients/mutations to perform selection analyses (Supplementary Figure 7). Significant signatures of purifying selection in pre-WGD truncating mutations in LOH in essential genes (but not after WGD) were observed in skin cutaneous melanoma (SKCM), which is characterized by high mutation rates and high LOH (Supplementary Figure 6). In other cancer types (liver hepatocellular carcinoma -LIHC-, colon adenocarcinoma -COAD-, uterine corpus endometrial carcinoma -UCEC-, bladder urothelial carcinoma -BLCA-, cervical squamous cell carcinoma and endocervical adenocarcinoma -CESC-, head and neck squamous cell carcinoma -HNSC-, esophageal carcinoma -ESCA-), although not significant, a similar trend was also observed.

Leveraging mutations in LOH to identify novel cancer genes

Positive selection is a feature of cancer, reflected as an enrichment of non-silent mutations in cancer genes^19,31,32, and thus the dNdS ratio can be exploited to identify new cancer genes. Next, we explored dNdS values specifically for known cancer genes. We observed highly significant positive selection of truncating mutations prior to WGD in regions of LOH, followed by reduced positive selection after duplication (Figure 4a). These data may reflect strong selection for inactivation of both copies of tumor suppressor genes through LOH and mutation, prior to WGD.

We hypothesized that by specifically looking into genomic segments of LOH we could identify tumor suppressor genes that would not be otherwise detected. Moreover, we could explore whether selection intensities varied depending on whether the mutation co-occurred with LOH or the wild-type allele remained intact.

We obtained dNdS values for nonsense mutations at the gene level and compared the selection intensities for early mutations in genomic regions of LOH (Figure 5a, x-axis) compared to mutations not occurring in segments of LOH (Figure 5a, y-axis). Notably, two genes, with very high selection coefficients in LOH but lower in non-LOH regions in LUSC were TP53 and PTEN (Figure 5a, left). These are well-established tumor suppressor genes that, concordant with the “two-hit hypothesis”³³, appear to require disruption of both alleles. ZNF750 (Zinc Finger Protein 750), is subject to high positive selection in regions of LOH in LUSC and would have remained undetected with standard procedures. Notably, ZNF750 has previously been described as a lineage-specific tumor suppressor gene in squamous cell carcinoma³⁴, consistent with its occurrence as an early event in 87% (7/8) of LUSC tumors. Other genes that were only identified as significant when looking at mutations in regions of LOH were NOTCH1 and SMAD4 in LUSC (Figure 5b and Supplementary Figure 3). On the other hand, CUL3 (Cullin-3), which plays an important role in the ubiquitin-proteasome system³⁵, showed high signatures of selection in genomic regions not subject to LOH, possibly reflecting haploinsufficient activity in cancer, and, conceivably, the requirement for an intact wild-type allele.

Finally, we extended the analyses to the remaining cancer subtypes in TCGA. A total of 33 potential tumor suppressor genes were identified by limiting the analyses to segments of LOH (Figure 5b). Some of these are included in the COSMIC and/or Network of Cancer Genes (NCG6.0) databases as well-characterized cancer genes (RB1, PTEN, SMAD4, BAP1, SETD2, NOTCH1, BRCA1, NF2 and CDK12). In other cases, we re-discovered genes that have been identified in other cancer types or pan-cancer studies, but not specifically in that cancer type, like SMAD4 in CESC, CDH1 in BLCA, ZNF750 in LUSC, ARHGAP35, TUBA3C and PTCH1 in HNSC, SMAD3 in COAD and ITGAV in READ. Finally, 42% of these were additional novel cancer genes may be acting as tumor suppressor genes not included in the current release of COSMIC and NCG6.0: WWC1, FPF12, NT5DC3, NCLN, KRTAP19-5, GRIK2, GLRA1, FAXDC2, FAM19A3, CRYGC, CLEC4E, BC02, ARPP21 and AC061992.1.

As an additional approach for the identification of driver genes under different scenarios (mutations in LOH, non-LOH, all) we also implemented MutSigCV^26,36. Supplementary Figure 8 shows those genes with q-values<0.1. We detected six potential cancer genes when looking at LOH mutations that were not significant when assessing all mutations together. These include the well-established cancer genes RB1, NCOR, CDNK2A, FBXW7 and TSC1, as well as PCDHGA3 in LUAD.

Data processing

We analyze 93 patients from the first cohort of patients of NSCLC (61 LUSC and 32 LUAD) obtained through the lung TRACERx (TRAcking Cancer Evolution through therapy (Rx)) project and thoroughly described in¹⁶.

Additionally, raw .bam files for LUAD (n=398) and LUSC (n=325) samples from the TCGA repository were downloaded and processed through the TRACERx pipeline. Briefly, we used BWA-MEM to align the reads to the reference genome (build hg19). We used Platypus for SNP calling on the germline, and Varscan2 and MuTect for somatic mutation calling. Functional annotation of genomic variants was performed using ANNOVAR⁴³. Purity, ploidy and copy number profiles of tumor cells were obtained with ASCAT⁴⁴, using the matching germline data. Mutations in regions of LOH were timed as early or late based on the mutation and major allele copy number. Following a conservative approach, we considered early mutations those with mutation copy number >=1.75 and major allele copy number >=1.75. Mutations were classified as late if mutation copy number <=1.25 and major allele copy number >=1.75. Clonal mutations that could not be timed were classified as “unknown”. Mutations were defined as mutations in LOH if the minor allele copy number was <0.25. The WGD status for each tumor was obtained using the genome doubling algorithm described in³.

For the identification of new drivers across other cancer types, we downloaded the MC3 somatic mutation tables⁴⁵ from whole-exome sequencing data across the 33 cancer types from TCGA. Further processing of the samples was performed as described above.

Simulations

We performed simulations to model the viability of a cancerous cell over time and illustrate the potential of increasing DNA copies via genome duplication in the mitigation of damaging effects of mutations using an adapted version of the cancer evolution model developed by McFarland et al.^14,15. Briefly, this model simulates cancer progression as a stochastic system of birth and death events using a standard Gillespie algorithm, and is defined by the mutation rate, the target sizes for drivers and passenger mutations, and the effect on the fitness of drivers and passenger mutations¹⁴. We have extended this model incorporating WGD events. The viability of the population of cells varies by a function of the mutations that the cell has accumulated over time, which exert deleterious or beneficial effects on cell fitness to varying degrees, depending on the type of genes where they are accumulated. Passenger alterations in regions of haploid LOH decrease the fitness by increasing the time to the next birth. Conversely, driver alterations decrease the time to the next birth.

For all the simulations we used the following fixed parameters based on the simulations by McFarland et al.¹⁴: starting population size (No) 1000 cells, 3000 driver loci (Td) with a maximum of 10 driver mutations in a cell, selection coefficient of the drivers (sd) 0.15. Different selection coefficients of passenger mutations (s_p) were tested, in the range of 0-0.01 (Table S1). We also explored different passenger alterations rates.

The effects of WGD events were incorporated by reducing the deleterious effects of passenger mutations by half. Note that given we are considering passenger mutations in regions of haploid LOH, we are likely underestimating the buffering effect of WGD.

The probability of WGD in a given cancer cell was defined by a normal distribution: if(rnorm(1, 0) > qnorm(0.001, lower = FALSE)). Every cell (and therefore its subsequent daughter cells) was only allowed to accumulate one WGD event. Additionally, we incorporated a fitness cost (range: 0-1) associated with surviving the WGD event that was maintained during fixed period.

In each simulation, the tumor progressed until the cell population size reached 20,000 cells - happening when a lineage reaches fixation in tumor progression very quickly - or until 2,000 generations in the cases where we don’t reach the maximum number of cells but the final population reaches an equilibrium in tumor size and composition. At that point we calculated the proportion of WGD cells in the tumor. An average of 16 iterations were run for each fitness cost.

Negative selection in simulations was evaluated by comparing simulations where there was no fitness cost associated with passenger alterations, i.e. s_p = 0, to those where s_p >0. For each s_p value, the median number of mutations at 1,000 generations was divided by the median number of mutations when s_p = 0. This ratio is equivalent to a dNdS ratio, such that a value <0 indicates positive selection, while a value <0 indicates a depletion of mutations and negative selection.

Copy number simulations

In order to explore whether there is less haploid LOH than would be expected by chance, we performed simulations based on copy number changes. For each of the WGD tumors in LUSC+LUAD from TCGA we determined the observed LOH and haploid LOH proportion and the total copy number gains and losses at the chromosome arm level. Using the total copy number gains and losses at the chromosome arm level as probabilities for gains and losses we simulated 10,000 random copy number scenarios and compared the simulated haploid LOH proportions with the observed ones. Essential genes scores per Chromosome arm were obtained from³⁰.

Selection tests

Selection across LUSC and LUAD was quantified using the dNdScv R-package¹⁹, an implementation of the traditional dNdS ratio, adapted and refined for cancer genomes. In all cases we used the trinucleotide substitution model with 192 rate parameters and default parameters, removing ultra-hypermutator samples (3,000 maximum number of coding mutations per sample) and limiting the analyses to 3 mutations per gene per sample. dNdS ratios were quantified for missense, nonsense and essential splice mutations, both globally and at gene level.

Global dNdS values for nonsense mutations were calculated for different lists of genes: a) all genes, b) cancer genes, which includes a compilation of LUSC and LUAD-specific tumor suppressor genes and oncogenes described in the literature^19,23–26 and c) essential genes, which includes the essential genes reported by Blomen et al.,²², identified by using extensive mutagenesis in haploid human cells (1,154 genes). Different lists of essential genes obtained through CRISPR-based systems were also evaluated for comparison^28,29. Essential genes from Hart et al.²⁸ included the 1,580 hits observed in three or more of the five cell lines. In the case of Wang et al., (2015), we included the genes with adjusted p-values <0.05 in three or more of the four cell lines used (a total of 1,282 genes). In all the analyses, cancer genes (present in the cancer gene census) were removed from the essential genes list.

Selection tests were performed on different subsets of mutations considering the timing (pre-WGD vs post-WGD mutations) and the LOH presence (LOH vs nonLOH) in WGD genomes.

Additionally, we calculated dNdS values for early mutations in LOH in essential genes inferred by NSCLC CRISPR analyses (Achilles) available at the DepMap portal https://depmap.org/portal/. Genes were ranked by the average gene-knockout effect (a measurement of the consequences of gene deletion on cell viability) across 283 samples annotated as lung cancer for the primary disease and took lists of top genes of different lengths ranging from 500 to 7,000.

Identification of driver genes

Maximum-likelihood dNdS estimates (MLE) at the gene level obtained with the dNdScv package were used to identify potential driver genes under positive selection. We compared the MLE values for nonsense mutations (wnon) in regions of LOH where the mutation was present on all copies of the remaining allele vs non-LOH in WGD and nWGD tumors combined. As an additional method to identify driver genes we used MutSigCV³⁶, comparing the adjusted p-values for mutations in LOH, non-LOH and all mutations combined. In order to reduce false positives and multiple testing correction we filtered the genes that are expressed at low levels (below the median expression value across all the genes and patients) in that particular cancer type. RNAseq-based normalized expression data for each cancer type was obtained from the TCGA data portal.

Validation of the timing approach using mutation and copy number data from diploid and tetraploid clones from the HCT-116 cell line

In order to demonstrate that we can accurately quantify which mutations and copy number losses occur pre and post the WGD event we used mutation and copy number data from 2 diploid and 4 tetraploid subclones derived from diploid human colon carcinoma HCT-116 cells, isolated at different passages (4 and 50 passages)³.

For each sample exome capture was performed on 1-2 ug DNA isolated from genomic libraries with median insert size of 190 bp, using a customized version of the Agilent Human All Exome V5 kit, according to the manufacturer’s protocol (Agilent). Samples were 100 bp paired-end multiplex sequenced on the Illumina HiSeq 2500 at the Advanced Sequencing Facility at the Francis Crick Institute. The data was aligned to the reference human genome (build hg19) and somatic mutations were then obtained with Varscan⁴⁶ and MuTect⁴⁷ as previously described¹⁶, using the parental diploid clone as the germline. Copy number profiles were inferred with PICNIC using genotype chip data, as previously described³. The early (pre-WGD)/late (post-WGD) classification of the mutations was assessed using the same procedure as above.

Statistics

R version 3.3.1 was used to analyze the data. No statistical tests were used to predetermine the same size. Tests involving comparison of distributions were done using ‘t.test’ using the unpaired option, unless otherwise stated. Confidence intervals for dNdS analysis was obtained using the dNdSCV package¹⁹.