Mutational load of the mitochondrial genome predicts pathological features and biochemical recurrence in prostate cancer

Prostate cancer patient clinical and pathological characteristics

All 115 patients presented with prostate cancer and elected for surgical removal of the prostate gland, with a mean of 107 months (range 24-150 months) follow-up. Patient and pathological characteristics may be found in Table 1. In brief, the mean age at radical prostatectomy was 61 years (range 46 to 75 years), with pre-operative PSA levels ranging from 2.43 to18.6 ng/ml. Since our cohort is biased towards more aggressive disease presentation, defined by a Gleason score ≥7 (89/115, 77%), a substantial number of patients (41/115, 36%) had relapsed within three years (range 1 to 126 months) post-surgery. Radical prostatectomy was non-curative in eight patients presenting with advanced stage prostate cancer, two with reported mortality, while eight patients have advanced to metastatic disease. Pathological reassessment of the core biopsy from which DNA was extracted for this study corresponded with the total and primary Gleason score at surgical diagnosis in 46% (53/115) and 71% (82/115) of cases respectively, further detailed in Supplementary Table 1, while tumor purity estimations averaged 57% (range 5% to 85%).

Mitochondrial haplogroup distribution and prostate cancer risk association

Previous studies have suggested that inherited mitochondrial genome variation which defines population-specific haplogroups, may be associated with prostate cancer risk. Specifically, a study of 221 white North American men showed increased risk of haplogroup U for developing PCa [odds ratio 1.95] [24]. Sequencing complete mitochondrial genomes from blood, we found the haplogroup frequencies did not deviate significantly (p = 0.95 by two-sided Fischer’s exact test) from that expected for an Australian European-ancestral population [25] with haplogroup H most commonly represented at 43% and 44% for expected and observed frequencies respectively, followed by U (both 14%), T (both 9%), J (11% versus 9%) and K (8% versus 7%). Haplogroups were also not associated with specimen Gleason score (p = 0.81 by two-sided Fischer’s exact test). Patient specific haplogroups are depicted in Supplementary Table 1. Our study is in agreement with a more recent report that suggests that common European-derived mtDNA haplogroups are not correlated with increased prostate cancer risk [26,27] while the potential association between the earliest diverging L0 African-derived mtDNA haplogroups and aggressive prostate cancer [22] remains to be confirmed.

Spectrum and frequency of acquired prostate cancer specific mitochondrial genome variation

Sequencing mtDNA from matched tumors, we identified a total of 76 somatic single nucleotide variants (SNVs) in 50 patients (43.5%) (Table 1). A single variant recurred in three patients at position mt.16093, resulting in 74 unique variants (see Supplementary Table 2 for full list of variants). While sequencing was performed using the Ion Torrent PGM instrument, we used the Illumina HiSeq X Ten platform to validated 11/12 (91.7%) unique SNVs within 16 randomly selected tumor-normal sample pairs. The single SNV not identified in the HiSeq data represented a loss of a heteroplasmic low allele frequency (0.17) SNV and is most likely to be a discrepancy in variant calling between the platforms. We retained the SNV for our analysis. Compared with publicly available data, 91.9% (68/74) of the observed unique variants are novel to prostate cancer, 86.5% (64/74) novel to any cancer, and 77.0% (57/74) novel to any disease. Although distributed across the entire mitochondrial genome, we found a predominance of variation within the control region and ribosomal RNA genes (Figure 1A). When correcting for length, we found a greater, though not significant, frequency of variation within the non-protein coding than protein coding regions (p = 0.013 by two-sided Fischer’s exact test) (Figure 1B). This proved to be significant for the D-loop and the tRNA regions (p = 0.0044 by two-sided Fischer’s exact test), as previously reported [16]. Of the 41 SNVs mapping to protein coding regions, 26 directly alter an amino acid (non-synonymous mutation) of which 21 have biological potential defined by termination of translation (n=2) or by PolyPhen-2 computational prediction as either possibly (n=2) or probably damaging (n=17). The gene coding for protein Cytochrome b was affected most by damaging variants, both in absolute numbers and when corrected for gene length. Somatic SNVs further showed a predominance of pyrimidine nucleotide transition events, 54% C>T and 24% T>C in line with what was reported previously [20]. To determine the potential of SNVs to impact mitochondrial function (>70% variant frequency), we corrected for tumor purity estimates, resulting variant frequencies ranged 0.17-1.79, (average 0.75). Following this adjusted variant frequency the cumulative variant frequency (CVF) in each patient with variants ranged 0.28-3.25 (average 1.42). Investigating a previously reported ~3.4 kb deletion variant, long range PCR did not show this deletion to be detectable in any of our patients [28].

Correlating somatic mtDNA variation with clinical presentation and outcomes

While we found no single somatic mtDNA mutations common between patients, we observed a diagnostic and prognostic correlation with the total somatic variant burden and for the CVF per patient. Specifically, an increase in total number of somatic mtDNA SNVs and CVF was significantly associated with high-risk disease defined as an increased Gleason score at radical prostatectomy p = 0.0132 (Figure 2A) and p = 0.011 (Figure 2B), respectively, and within the sampled biopsy core p = 0.0013, (Figure 2C) and p = 0.0143 (Figure 2D), respectively, with p-values estimated from ANOVA models. Analyses for subgroups of variants showed a significant predictive value for number of non-synonymous variants (p = 0.034 by ANOVA) and synonymous and non-coding variants (p = 0.018 by ANOVA). We point out that, when adjusting for multiple testing using the conservative Bonferroni correction, all six tests drop below the significance threshold at the familywise alpha level of 5% (α = 0.05/7 = 0.0071). However, we note these tests are not independent, with two pairs of highly correlated variables, namely: number of SNV and CVF, and Gleason Scores at prostatectomy and of the biopsy cores.

When investigating the correlation between Gleason score and tumor purity we found a strongly significant correlation (p = 1.6e-5 by ANOVA for the linear model of: Gleason Score ~ Purity), consequently, the number of SNVs also correlates significantly with tumor purity (p = 0.003 by ANOVA for linear model of: SNV number ~ Purity). Additionally, patients with one or more somatic SNVs were significantly more likely than patients without any somatic mtDNA variants to present with a disease relapse (p = 0.012, hazard ratio (HR) 2.17 (confidence interval (ci) 1.17 – 4.04) by log rank test, Figure 3A), with risk for relapse increasing when considering 2 or more SNVs (p < 0.0001, HR 3.82 (ci 1.54 – 9.46) by log rank test (Figure 3B). As Gleason score is currently the most significant predictor for prostate cancer, we determined the power of total mtDNA variation burden to predict disease relapse in our study either alone and in combination with Gleason score, suggesting that number of mtDNA SNVs slightly improves relapse prediction in combination with Gleason score (p = 0.085 by DeLong test for area under curve (AUC)) (Figure 4).

Patient selection and sample preparation

Patients were recruited and consented and biospecimens biobanked according to St Vincent’s Hospital Human Research Ethics approval #HREC/12/SVH/323 at time of surgery. Study inclusion criteria was as follows; (i) a clinicopathological confirmed diagnosis of prostate cancer, (ii) no pre-operative treatment received (including hormone deprivation therapy, chemo or radio-therapy), (iii) a follow up of more than 24 months post radical prostatectomy, and (iv) the availability of a fresh frozen biopsy core with a matching blood sample. Biopsy cores were processed using cryostat-assisted frozen section sample preparation for pathology review (minimally two pathologists JG, RG, and/or JGK) to ensure tumor presence, evaluate the core pathology for percentage adenocarcinoma cells and Gleason score, and maximize for tumor purity during DNA extraction. DNA was extracted from matching blood and prostate tissue using standard methods (Qiagen) and made available under additional site-specific St Vincent’s Hospital Human Research Ethics approval #HREC/15/SVH/227 for whole mitochondrial genome sequencing. Based on these criteria a total of 115 patients were included in this study with clinic-pathological data related to age, PSA level and Gleason score at diagnosis, as well as follow-up with clinical and biochemical relapse data summarized in Table 1 and further elaborated in Supplementary Table 1.

Mitochondrial genome sequencing

To ensure mitochondrial genome specificity and thereby avoiding amplification of nuclear copies of mitochondrial-derived DNA (NuMTs), we performed a dual-amplicon long-range pre-amplification of the mitochondrial genome prior to sequencing on the Ion Torrent PGM instrument, as previously described [22]. In brief, after quantification, the 7.2 kilobase (kb) and 9.7 kb amplified products were combined in a 7:13 ratio, 200 bp libraries prepared using the Ion Xpress Plus Fragment Kit and Ion Xpress™ Barcode Adaptors (ThermoFisher), and six to eight matched normal-tumor pairs pooled and run on a 318v2 Chip. We generated per patient an average of 80 Megabytes (MB) of data, with an average coverage depth of more than 2000 X per mitochondrial genome.

Variant calling and annotation

Sequencing data were analyzed using the Ion Torrent suite v5.0.2.1. Specifically, sequencing reads were quality trimmed, aligned to the revised Cambridge Reference Sequence (rCRS, accession number: NC_012920) [32] and variants called using the Torrent Variant Caller v5.0.2.1 with customized parameters optimized for the mitochondrial genome. Somatic variants were identified by patient matched tumor-normal (blood) mitochondrial genome comparison, using an allele frequency difference greater than 10%. Variant annotation was performed using ANNOVAR [33] and functional impact was predicted using PolyPhen-2 [34]. In addition to single nucleotide variants (SNVs) and indels (small insertions and deletions), we used a previously published long-range 3.8 kb targeted amplification method [22] to screen for the 3379 bp mtDNA deletion suggested to be associated with prostate cancer in a north American cohort [28].

Variant validation

A total of 16 normal-tumor pairs (14% of the study) were randomly selected to undergo whole genome paired-end 2 x 150 bp sequencing using the Illumina HiSeq X Ten platform. A mean coverage of 44X and 71X was generated for the blood and tumor, respectively, resulting in mitochondrial genome coverage of roughly 3000 X and 14000 X, respectively. Mitochondrial sequences were trimmed and filtered and aligned to rCRS using bwa-mem v0.7.12 [35]. SNVs and indels were called using Strelka [36] and manually compared to sample matched SNVs/indels using the Ion Torrent PGM. Although several somatic indels were called using the Ion Torrent data, none were validated and no additional indels observed using the HiSeq X ten data. Conversely, 11/12 Ion Torrent called SNVs were HiSeq X ten validated and therefore only SNVs were used in subsequent mtDNA somatic mutation analyses.

Clinical correlation and statistical analysis

Clinically relevant data was extracted using CanSto ®, an in-house clinical database manager. Presentation at diagnosis was defined by three Gleason score categories (< 7, equal to 7, and > 7). Disease relapse was defined in this study as; non-curative surgery or a rise in PSA >0.2 ng/mL in two consecutive measurements (biochemical recurrence) either with or without the presence of distal visceral or bone metastasis, and/or associated death. Fisher’s exact test was used to test for difference in expected versus observed haplogroups, for difference in presentation at diagnosis and for testing expected versus observed variants in mtDNA regions. ANOVA linear regression models were used to test for clinical correlation, either Gleason score at radical prostatectomy or for the specific biopsy core, with mtDNA somatic mutation rate defined as SNV number or CVF. When more than two statistical tests were performed on subgroups, Bonferroni multiple comparison correction was performed using an alpha of α = 0.05. Statistical analysis and figures were generated in R (R Foundation, Vienna, Austria) using the ggplot2 package [37], with receiver operating characteristic (ROC) curves and area under the curve (AUC) displayed and analyzed using pROC [38]. Prism 7.0 was used to draw and test survival differences using Kaplan Meier Curves.