Cancer mutational burden is shaped by G4 DNA, replication stress and mitochondrial dysfunction

Abstract

A hallmark of cancer is genomic instability, which can enable cancer cells to evade therapeutic strategies. Here we employed a computational approach to uncover mechanisms underlying cancer mutational burden by focusing upon relationships between 1) translocation breakpoints and the thousands of G4 DNA-forming sequences within retrotransposons impacting transcription and exemplifying probable non-B DNA structures and 2) transcriptome profiling and cancer mutations. We determined the location and number of G4 DNA-forming sequences in the Genome Reference Consortium Human Build 38 and found a total of 358,605 covering ∼13.4 million bases. By analyzing >97,000 unique translocation breakpoints from the Catalogue Of Somatic Mutations In Cancer (COSMIC), we found that breakpoints are overrepresented at G4 DNA-forming sequences within hominid-specific SVA retrotransposons, and generally occur in tumors with mutations in tumor suppressor genes, such as TP53. Furthermore, correlation analyses between mRNA levels and exome mutational loads from The Cancer Genome Atlas (TCGA) encompassing >450,000 gene-mutation regressions revealed strong positive and negative associations, which depended upon tissue of origin. The strongest positive correlations originated from genes not listed as cancer genes in COSMIC; yet, these show strong predictive power for survival in most tumor types by Kaplan-Meier estimation. Thus, correlation analyses of DNA structure and gene expression with mutation loads complement and extend more traditional approaches to elucidate processes shaping genomic instability in cancer. The combined results point to G4 DNA, activation of cell cycle/DNA repair pathways, and mitochondrial dysfunction as three major factors driving the accumulation of somatic mutations in cancer cells.

Keywords: Cancer mutations; G-quadruplexes; Genome instability; Mitochondrial dysfunction; Replication stress; Translocation breakpoints.

Figure 1.. G4 DNA structures are readily detected in cell nuclei.

Confocal microscopy of 293T, HAP1 and Hela cells stained with DAPI (blue) for nuclear DNA, a G4 DNA-structure specific antibody (red) and with Phalloidin for cytoplasmic cytoskeleton (green) display nuclear colocalization of G4 DNA structural foci with chromosomal DNA.

Figure 2.. G4 DNA elicits translocations in cancer genomes.

Panel A, plot of fraction of translocation breakpoints or control genomic positions (y-axis) occurring within 5 bp of hg38 genomic coordinates mapping G4 DNA-forming sequences (x-axis). Brown, translocation breakpoints in cancer genomes. Green, 0.2, 0.4, 0.6, 0.8, 1.0 and 10.0 kb before (left) and after (right) cancer translocation breakpoints. Gray, random genomic coordinates. Panel B, number of translocation breakpoints (mean ± SD) plotted for patients with (brown) and without (green) breakpoints at G4 DNA-forming sequences. Panels C and D, bar graphs of fraction of patients without (panel C) or with (panel D) translocation breakpoints at G4 DNA-forming sequences harboring pathologic mutations in the top 20 cancer-mutated genes. Panel E, net fraction of cancer patients with translocation breakpoints at G4 DNA-forming sequences harboring pathologic mutations in the top 10 most mutated genes; only hits with a Pα_0.05 > 0.8 were recorded. Panel F, bar graph of percent transposable elements (TE) harboring G4 DNA-forming sequences in hg38 (orange) and at cancer translocation breakpoints (blue); all seq, all G4 DNA-forming sequences at SVA elements; main seq, most common G4 DNA-forming sequence at SVA retrotransposons (see Panel G). Panel G, list of most common G4 DNA-forming sequence at SVA elements (top); COSMIC ID, COSMIC tumor identification number; Tumor type, PRC, prostate cancer; OVC, ovarian cancer; UTC, uterine cancer; PAC, pancreatic cancer; Hg38 coor, genomic coordinate of translocation breakpoint within G4 DNA-forming sequence; TE, SVA lineage; Tot trans BP, total number of translocation breakpoints in the tumor sample.

Figure 3.. Gene expression profiles correlate with cancer somatic mutations.

Panel A, Plot of regression coefficients (R, x axis) vs. P-values (y axis) for correlations between gene expression (all genes) and somatic mutations for patients with CHOL (black) or BRCA (red). Panel B, S-plots of P-values for the correlations between gene expression (all genes) and somatic mutations for 32 TCGA datasets. Panel C, list of top 10 genes whose mRNA levels (expression) were most strongly correlated with somatic mutations. GO term, selected gene ontology term from the human gene database GeneCards (https://www.genecards.org); Corr, P-value, P-values from panel B; KM P-value, P-value for the Kaplan-Meier plot; COSMIC CGC, whether or not the gene is listed as cancer-promoting gene in the COSMIC cancer gene census. Panel D, Kaplan-Meier survival curves for LUAD patients with low (red) gene expression for PIGR and SPATA18 relative to LUAD patients (blue) with the other 3 combinations (high PIGR high SPATA18; low PIGR high SPATA18; high PIGR low SPATA18). Panel E, box plot of mRNA levels (y axis) in the tumor (blue) and normal (red) samples for the 10 genes with the strongest positive correlations between gene expression and somatic mutations for patients with LUAD. ***, P-value <2 × 10⁻¹⁶. Panel F, number of tumor types (y axis) in which gene expression for a given gene (x axis) was higher in the tumor than in matched normal control tissues (green, 15 total tumor types tested) and the number of instances (orange) in which the P-value for the difference was <2 × 10⁻¹⁶.

Figure 4.. Altered gene expression in key pathways decreases patient survival.

Panel A, box plot of the correlation coefficient R for the coexpression of MYBL2 vs. CENPA, KIF2C or KIFC1 in patients for each of the 32 TCGA tumor types. Panels B and C, line plots for the normalized Rsem gene expression values for MYBL2 vs. KIFC1 in LGG patients (panel B) and vs. KIF2C in BRCA patients (panel C). R, regression coefficient; P, P-value. Panel D, list of tumor types with worse survival for high MYBL2 gene expression levels and P-values for the respective Kaplan-Meier survival curves. Panel E, Kaplan-Meier survival curves for MESO patients with low (red) and high (blue) gene expression levels for MYBL2. Panel F, list of tumor types in which expression of the succinate dehydrogenase complex genes was lower in the tumor than in the matched control tissues and the corresponding P-values. Panel G, Kaplan-Meier survival curves for KIRC patients with low (red) and high (blue) SDHD gene expression.

Figure 5.. Mutations are linked to defects in DNA repair gene expression.

Panels A - C, bar graphs of P-values adjusted for multiple testing (Benjamini correction) for genes enriched in KEGG pathways from the pool of top 270 genes (for each tumor type) that displayed positive correlation between mRNA levels (expression) and cancer somatic mutations. Panels A and B group tumors that displayed common terms, whereas panel C shows singletons of tumor-specific KEGG pathway enrichment. Panel D, graph of regression coefficients (y axis) for DNA repair genes (x axis) displaying P-values <1 × 10⁻⁵ for the correlations between gene expression and cancer mutation loads in patients with LUAD; shaded areas highlight genes in the same pathway; BER, base excision repair; MMR, mismatch repair; NER, nucleotide excision repair; HR, homologous recombination; FA, Fanconi anemia.