Kataegis co-localises with MYC translocation breakpoints

We were also able to detect a mutational signature, kataegis, where regional clustering of mutations can be indicative of somatic genomic rearrangements¹⁴. It is difficult to detect kataegis using exome sequencing but the tiled regions surrounding the IGH, IGK, IGL and MYC loci could be used to detect it. By creating rainfall plots we were able to discern samples with regional hypermutation. As expected, the Ig loci contained clusters of hypermutation, but these were not enriched for C>T or C>G mutations within a TpCpH trinucleotide context and as such are not caused by kataegis. We found the hallmarks of kataegis in 15 samples (3.2%), where there was enrichment for TpCpH mutations with an inter-mutational distance <1 kb (Figure 2). Of these 15, 9 were found in the tiled region surrounding MYC and others were detected on chromosomes 1, 10, 11, 16, and 17. Kataegis was co-localised with copy number abnormalities in 12 of the samples. Two of the samples with kataegis surrounding MYC also had an inter-chromosomal translocation at MYC involving either IGK or IGL. Interestingly, the partner chromosomes also showed signs of kataegis e.g. in the t(2;8) kataegis was found at IGK and MYC and in the t(8;22) kataegis was found at MYC and IGL (Figure 2). The pattern of mutations clustered around the translocation breakpoint and according to the cancer clonal fraction (CCF) were present in all cells, indicating that kataegis most likely occurs at the same time as the translocation. APOBECs are thought to be involved in the generation of kataegis and as such this co-localisation is indicative of APOBEC involvement in the generation of MYC breakpoints.

Open in a separate window

Figure 2

Kataegis in myeloma

Kataegis at 8q24 coincides with translocation breakpoints and occurs on both chromosomes where a translocation has occurred. Examples of two samples are shown. For each sample the top left plot shows the distances between mutations and are colored on a chromosomal basis according to UCSC coloring. The top right panel show the same data but is colored by mutation type, as per the key. The bottom panels show the location on chromosome 8 and the partner translocation chromosome (22 or 2, respectively) where kataegis is found with genes or Ig loci segments indicated in green or cyan, respectively. The arrows indicated the position of the translocation breakpoint.

Hypermutation of Translocation Partner Oncogenes

In agreement with previous studies¹⁵ we find that CCND1 is mutated, and this was seen in 10 patients, all of whom have a t(11;14) (Supplementary Table 3). All mutations occurred in the first exon of CCND1 and all but one were located outside of the cyclin box fold domains. Five patients had multiple mutations and the K50R mutation was detected in three samples, the K46N mutation was seen in two samples but all of the other mutations we describe were unique. There was no association of mutation in CCND1 with type of breakpoint (RAG or switch mediated) or allele of the variant associated with the t(11;14). There was a significant association with mutation of CCND1 and the distance to the translocation breakpoint, where samples with a breakpoint closer to CCND1 were more likely to have mutations within CCND1 (p=0.03; Supplementary Fig. 3). There was also an association of mutated CCND1 and a poor prognosis when compared with non-mutated t(11;14) patients (Overall survival median of 20.2 months versus not reached, log rank p=0.005; Figure 3C) in the Myeloma XI trial. This effect was also seen when comparing patients with mutated CCND1 to all other patients in the cohort and was an independent prognostic factor for both PFS and OS. To determine the significance of this result we sequenced the first exon of CCND1 in 102 t(11;14) samples from the Myeloma IX trial and found mutations in a further 10 samples (9.8%; Supplementary Fig. 4). There was no effect on survival in the Myeloma IX dataset with mutations in CCND1. We also examined the allelic frequency of the variant associated with the t(11;14), rs9344, in the Myeloma XI dataset and found that in agreement with our previous observations⁸ the G allele is significantly associated with the translocation (Supplementary Table 2).

Open in a separate window

Figure 3

Mutations in translocation partner oncogenes

(A) Non-synonymous mutations in translocation partner oncogenes are depicted along with the translocation group they occur in. (B) The cancer cell fraction (CCF) in which the mutations occur. (C) Mutation of CCND1 in t(11;14) samples results in decreased overall survival in Myeloma XI samples. The box and whisker plot shows the 10 and 90 percentiles (whiskers) and the median and 25 and 75 percentiles (box).

Given the association of mutated CCND1 in the t(11;14) we examined the other translocation groups and their associated partner chromosome oncogenes. We found that FGFR3, MAF and MAFB all had mutations which were restricted to the relevant translocation group (16.9% mutated FGFR3 in t(4;14); 12.5% mutated MAF in t(14;16); and 25% mutated MAFB in t(14;20); Figure 3A). There were no mutations in CCND3 (in 6 t(6;14) patients) and the mutations in MMSET were not in t(4;14) patients. In contrast to the poor prognosis associated with mutations of CCND1 in the t(11;14), no poor prognosis was associated with mutations in FGFR3, MAF or MAFB. However, the latter two sample sets may be too small to address this question adequately.

The mutated oncogenes were all on the der(14) and most likely reflect somatic hypermutation events mediated by AID, an member of the APOBEC family, which would normally affect the V(D)J rearrangement upstream of the IGH constant regions. After the translocation event, MMSET is on the der(4) and is, therefore, unlikely to be mutated by this mechanism. To determine that AID may be responsible we examined the transition/transversion ratio and the proximity of the mutations to the AID consensus motif, WRC. In agreement with previous studies¹⁶ the mutations in the IGH translocation partner genes are 58.8% in A:T bases, of which 50% are transversions (Supplementary Fig. 5), and in close proximity to WRC motifs (Supplementary Table 3) indicating that AID is involved in this mutational process. As AID is associated with hypermutation of the V(D)J and switch regions of the highly expressed IGH locus we examined the expression of the translocation partner oncogene between those with or without a mutation, but found no significant differences (Supplementary Fig. 6).

The CCF in which the mutation was found differs by translocation group (Figure 3B). In the t(11;14) the mutations were founder events, present in all of the cells whereas in the t(14;16) and t(14;20) the mutations were in only ~50% of the cells indicating they are obtained later than the translocation themselves. In the t(4;14) the mutations in FGFR3 can be either clonal or subclonal indicating that these mutations can develop at the same time as the translocation or at a later time point.

The mutation patterns for each of the der(14) oncogenes differ. In the t(11;14) mutations in CCND1 occur solely near the N-terminus of the protein. In MAF and MAFB the mutations are constrained in the basic-leucine zipper domain at the C-terminal of the protein. The focussed mutation profile seen in CCND1, MAF and MAFB are indicative of activating mutations. The mutations in FGFR3 are dispersed through several domains but have also been described as mutations that can activate the RAS/MAPK pathway in urothelial cancer and myeloma (Supplementary Fig. 7) ^{17, 18}.

APOBEC Mutations are enriched in ‘maf’ translocation groups

It has previously been shown that mutations can be described in a specific trinucleotide context and that in myeloma there are two signatures that predominate¹⁰. We are able to recapitulate these two different signatures, which consist of a generic signature comprised of enrichment of C>T transitions in a CpG context, Signature A (Figure 4A) and a second signature, Signature B, in which there is enrichment for C>G and C>T, especially in a TpCpA context. Signature B is hypothesised to result from aberrant APOBEC activity, where the APOBECs enzymatically modify single-stranded DNA. AID is a member of the APOBEC family and is involved in class-switch recombination and somatic hypermutation in B cells.

Open in a separate window

Figure 4

Analysis of mutation context identifies two signatures in myeloma

(A) Mutation context identifies two signatures in myeloma, Signature A and Signature B. (B) t(14;16) and t(14;20) samples have significantly more mutations than other translocation groups. (C) The mutational context split by translocation group identifies t(14;16) and t(14;20) with more mutations which make up Signature B. (D) Stacked bar chart showing the percentage contribution of the two signatures identified by NMF in each sample, ordered by translocation group. The box and whisker plot shows the 10 and 90 percentiles (whiskers) and the median and 25 and 75 percentiles (box). P-values determined by ANOVA.

We noted that the t(14;16) and t(14;20) have a statistically significant higher number of mutations per sample compared to the other translocation sub-groups (2-sided t-test p=1.65×10⁻⁵), Figure 4B, and that the Signature B (APOBEC) related context of mutations in the t(14;16) and t(14;20) was significantly higher than other translocation groups (T(C>T)A, p=9.1×10⁻⁶; T(C>T)T, p=0.0014; T(C>G)A, p=0.001; T(C>G)T, p=0.0064), Figure 4C. Collectively, mutations in these trinucleotide contexts comprise a mean of 28.7% and 21.3% of mutations in t(14;16) and t(14;20) samples, respectively (compared to t(4;14) 6.5%, t(6;14) 6.2%, t(11;14) 5.8%). We examined the proportion of each Signature present in the translocation sub-groups and found that there is a significant enrichment for Signature B (APOBEC) mutations in the t(14;16) and t(14;20) (0.56, p=2×10⁻¹⁶; 0.44, p=8.26×10⁻¹¹, respectively) compared to the t(4;14), t(6;14), t(11;14) and hyperdiploid samples (0.094, 0.096, 0.074, 0.098 and 0.078, respectively) (Figure 4D). These data indicate that the ‘maf’ translocation groups are largely characterised by APOBEC signature mutations and have a higher mutation load than the other translocation groups.

In order to determine if there are some samples in the other translocation groups which also have an APOBEC signature we assigned each sample to either Signature A or Signature B depending on the proportion of mutation type in each sample. This generated clusters of samples whose mutations are either mostly Signature A (cluster A) or Signature B (cluster B), Figure 5A. Here we find that the t(14;16) and t(14;20) cases comprise 66.6% of cluster B but only 1.3% of cluster A. In line with the proportion of ‘maf’ samples in Cluster B the number of mutations in this cluster is significantly higher than in Cluster A (mean 295.44 vs. 127.22, two-sided t-test p=1.18×10⁻¹⁵; Figure 5B). Cluster A is comprised of 445 patients and cluster B 18 patients, indicating that Signature B only affects 3.9% of patients. However, when we performed survival analysis of these patient clusters we find that there is a significant effect on overall survival (log rank p=0.02; Figure 5C). Due to the low number of samples with ‘maf’ translocations and the interaction of the translocation, APOBEC signature and mutational load this effect on survival is not significant in a multivariate analysis as it is not possible to delineate whether this effect on survival is due to any single one of these markers alone and it remains more likely that the impact of three abnormalities is linked.

Open in a separate window

Figure 5

Myeloma mutations can be categorised as belonging to Signature A or Signature B

(A) Samples which mostly have Signature B mutations dominate Cluster B and are enriched for t(14;16) and t(14;20). P-value determined by a linear model. (B) Samples in Cluster B have more mutations than those in Cluster A, p-value determined by ANOVA. Patients in Cluster B have a significantly worse overall survival, p-value by log rank test (C). The box and whisker plot shows the 10 and 90 percentiles (whiskers) and the median and 25 and 75 percentiles (box).

Multiplex Ligation-Dependent Amplification

50 ng of sample DNA was subjected to MLPA reactions using the SALSA P425-B1 probemix developed by MRC-Holland (Amsterdam, The Netherlands).³¹ The probemix contained 46 probes for copy number analysis of the major prognostic markers in myeloma including: 1p32-31 (FAF1, CDKN2C, PPAP2B and DAB1), 1p21, 1p12 (FAM46C), 1q21.3 (CKS1B), 1q23.3, 5q31.3, 12p13.31 (CD27, CHD4), 13q14 (RB1, DLEU1, DLEU2, DIS3), 16q12 (CYLD), 16q23 (WWOX), and 17p13 (TP53). In addition, 11 reference probes were also included, analyzing various autosomal chromosome loci which are only infrequently involved in CNA in myeloma. A second panel was used to interrogate the MYC locus, which contained probes for exons 1 and 3. The MLPA reaction, including internal quality controls and negative controls, was performed according to manufacturer’s instructions. The PCR-products were analyzed using an ABI 3730 DNA analyzer (Life Technologies, Paisley, UK) and Coffalyser software (MRC-Holland, Netherlands). After intra-sample and inter-sample normalizations, copy number at each locus was estimated.³² In summary, values ≥1.2; between 1.19 and 0.71; and ≤0.7 were considered as gain, normal, and deletion, respectively. If multiple probes per gene were examined a consistent result between them was required (2 probes with same result or majority of probes with same result depending on locus/chromosome examined).³³

Gene Expression

Translocations were defined using the TC classification based on expression of the Ig-partner oncogene. In addition, MYC expression was assayed (Hs00153408; Life Technologies, Paisley, UK) using qRT-PCR and relative quantitative (RQ) values were compared with the endogenous control (GAPDH).

Exome Sequencing and Analysis

DNA from both tumor and peripheral blood samples were used in the exome capture protocol as previously described.³⁴ Briefly, 200 ng DNA was fragmented on a Covaris E-Series using the settings suggested in the Sureselect^XT Target Enrichment System for Illumina Paired-End Sequencing Library protocol (version 1.5). DNA was end-repaired, A-tailed and adaptors ligated using the NEBNext DNA library prep master mix set for Illumina (New England Biolabs, Hitchin, UK). Adaptor-ligated DNA was amplified using the NEBNext High-fidelity PCR master mix (New England Biolabs) using 7 PCR cycles and 750 ng DNA hybridized against RNA baits overnight (SureSelect Human All Exon V5, Agilent Technologies). The RNA baits were designed against the human exome with the addition of custom baits that tiled the IGH, IGK, IGL and MYC loci in order to detect the major translocations in myeloma. After hybridisation the captured DNA was indexed and amplified using Herculase II fusion DNA polymerase (Agilent Technologies) for 8 PCR cycles. Four exome samples were pooled and run on one lane of a HiSeq 2000 (Illumina, Hinxton, UK) using 76-bp paired end reads.

Data Quality Metrics and Processing

FastQC (v0.10.0) was used for basic quality control of Illumina paired-end sequencing data. These files were aligned to the reference genome (GRCh37), using BWA (v.0.5.9) followed by Stampy (v.1.0.20) to improve gapped alignment. BAM files were recalibrated using GATK (v2.3.9) and deduplicated using Picard (v.1.85). Tumour and normal samples were realigned as pairs using the GATK indel realigner to improve indel call rates. Calls from a panel of 6668 highly variable SNPs within the exome capture were compared between the tumour and normal samples to confirm that they were correctly paired.

Somatic Mutation Calling

Single nucleotide variants (SNVs) were called using MuTect (v1.1.4). They were further filtered using the following criteria: minimum of 10x depth at that site in both tumour and normal BAMs, a minimum of 1 non-reference base call in both directions, a mean Phred quality score of greater than 26 for that base in the tumour sample, a mean mapping quality score of equal to or greater than 50 across all reads at that site in the tumour sample, and the site must be uniquely alignable according to the CRG alignablity tracks available from the UCSC genome browser and created using the GEM mappability tool (Fast Computation and Applications of Genome Mappability by Derrien et al, 2012). C>A|G>T SNVs likely to be oxidation artefacts created during library preparation were removed.³⁵

Short indels were called using the GATK Indelocator and furthered filtered according to the following criteria: a minimum of 1 non-reference base call in both directions, a mean Phred quality score of greater than 26 for that base in the tumour sample, a mean mapping quality score of equal to or greater than 50 across all reads at that site in the tumour sample, and the Site must be uniquely alignable. No more than 2 reads covering that site in the normal sample could contain any indel. A window of 21 base pairs centred on the first base of the indel was taken and had to conform to the following rules which remove low-complexity sequences: no homopolymers of greater than 6 base pairs and no dinucleotide may occur more than 5 times.

Deletions of whole exons (windows defined as the regions used in the Agilent exome capture) were found by comparing the read depth between the tumour and normal samples. The mean depth across the window was required to be greater than 0.2 of the median depth in the normal sample and less than 0.06 in the tumour sample, with the normal value being at least 8x greater than the tumour value.

Copy number across the exome was determined using Control-FREEC³⁶ utilising 500 base pair bins, each overlapping with the subsequent and previous 250 bp. A minimum average read depth of 50 was required in the control samples, with at least two neighbouring bins required to show copy number aberration to call a region as gained or lost.

All somatic events were annotated using both SnpEff (v3.1) and Oncotator (v0.4.2.2) with SnpEff providing the most deleterious interpretation regardless of transcript and Oncotator annotating only the canonical transcript. Significantly mutated genes were detected by providing all SNV and short indels to the MutSigCV (v1.4) algorithm. A q-value cut-off of 0.1 was used.

The proportion of tumour cells containing an SNV was calculated using the following equation: $\mathit{CCF}=\min (1,\frac{\mathit{CN}.r}{R})$

CCF = cancer cell fraction (proportion of cells containing the mutation)

CN = copy number at that site

r = number of reads containing the mutation at that site

R = total number of reads at that site

Translocations were called using Delly³⁷ and manually curated in the Integrative Genomics Viewer (IGV). Translocations were considered real when there were at least 10 supporting reads and if no translocations were also found in the peripheral blood sample. Simultaneously, translocations were determined by qRT-PCR using defined cut-offs for expression of each partner oncogene.³⁸ Results from both techniques were compared and a consensus translocation call generated. There was concordance between assays in 93.4% of cases. When assays were not concordant the data were examined and a decision made.

Nonnegative Matrix Factorization (NMF)

Nonnegative matrix factorization (NMF) using the NMF package in R was used to factorize a matrix of frequency of trinucleotide mutation contexts per sample, in order to identify underlying mutation signatures among all mutations in the 463 samples.^{10, 14} The algorithm was run for between two and 18 underlying signatures, with 10 runs at each number. Two signatures were selected as the optimum number of signatures based around a number of quality control metrics. Samples were clustered using the consensus hierarchical clustering procedure implemented in the NMF package.

Correlation studies

To determine correlation between cytogenetic abnormalities a probabilistic approach based on Bayesian interference was determined in R ³⁹ using the program "JAGS"⁴⁰ and its the R interface Bayesmed ⁴¹ as previously described ⁴². The probability of the observed data under the null hypothesis versus the alternative hypothesis or Bayes factor (BF) was computed. BF>1 was considered significant.

Survival analysis

Survival analysis was performed in R³⁹ using package survival^{43, 44} and coin^{45, 46} Differences between curves were tested for statistical significance using the log-rank test, with p<0.05 taken as the level of significance.

Cell lines and Cell culture

Sachi is a cell line derived from a t(14;20) primary plasma cell leukemia with increased levels of MAFB transcripts and was derived by Hanamura and colleagues (2001).^{47, 48} RPMI8226 is a cell line obtained from the ATCC with a t(14;16) translocation and increase levels of MAF transcripts.⁴⁹ RPMI8226 and Sachi were cultured in RPMI1640 (Invitrogen, Carlsbad, CA) containing 10% heat-inactivated FBS, and 4mM L-glutamine, as described previously.⁵⁰ 293T cells were cultured in Dulbecco Modified Eagle Medium (DMEM) containing 10% heat-inactivated FBS, penicillin (100 U/mL), streptomycin (100 mg/mL), 4 mM L-glutamine, and 1 mM sodium pyruvate. The Cells were maintained at 37°C and humidified with 95% air and 5% CO₂ for cell culture.

Silencing c-maf and mafB expression by short hairpin RNA

The knockdown of c-Maf and MafB genes were performed in myeloma cells.⁵¹ Briefly, the sequence for shRNA specific to human c-MAF gene ( 5′-CCGGTGGAAGACTACTACTGGATGACTCGAGTCATCCAGTAGTAGTCTTCCATTTTT-3′) and to MafB (5′-CCGGGCCTTGTCTTATGGTCAAATTCTCGAGAATTTGACCATAAGACAAGGCTTTTT-3′ ) in lentiviral vector were transfected into 293T cells and supernatant from culture media were harvested and concentrated. A control oligonucleotide sequence not matching any sequence in the human genome (5′-GATCCCCGACACGCGACTTGTACCACTTCAAGAGAGTGGTACAAGTCGGTCGTCTTTTTA −3′) was used as a control shRNA sequence (designated as shCont). The myeloma cells were infected with lentivirus supernatants for 48 hours and total RNA was isolated using TRIzol (Life Technologies) according to the manufacturer’s instructions.

The authors would like to thank all the patients and staff at centres throughout the UK whose participation made this study possible. The authors are grateful to all principle investigators for their dedication and commitment to recruiting patients to the study. The principal investigators at the four top recruiting centres were Dr Don Milligan (Heart of England NHS Foundation Trust), Dr Jindriska Lindsay (Kent and Canterbury Hospital), Dr Nigel Russell (Nottingham University Hospital) and Dr Clare Chapman (Leicester Royal Infirmary). The support of the Clinical Trials Research Unit at The University of Leeds was essential to the successful running of the study and the authors would like to thank all the staff including Helen Howard, Corrine Collett, Jacqueline Ouzman, Ana Quartilho and Alex Szubert. We also acknowledge The Institute of Cancer Research Tumour Profiling Unit for their support and technical expertise in this study.

This work was supported by a Myeloma UK program grant, Cancer Research UK CTAAC sample collection grants (A12136 and A17761) and a Cancer Research UK Biomarkers and Imaging Discovery and Development grant (A14261) as well as funds from the National Institute of Health Biomedical Research Centre at the Royal Marsden Hospital.