Mutational mechanisms shaping the coding and noncoding genome of germinal center derived B-cell lymphomas

Abstract

B cells have the unique property to somatically alter their immunoglobulin (IG) genes by V(D)J recombination, somatic hypermutation (SHM) and class-switch recombination (CSR). Aberrant targeting of these mechanisms is implicated in lymphomagenesis, but the mutational processes are poorly understood. By performing whole genome and transcriptome sequencing of 181 germinal center derived B-cell lymphomas (gcBCL) we identified distinct mutational signatures linked to SHM and CSR. We show that not only SHM, but presumably also CSR causes off-target mutations in non-IG genes. Kataegis clusters with high mutational density mainly affected early replicating regions and were enriched for SHM- and CSR-mediated off-target mutations. Moreover, they often co-occurred in loci physically interacting in the nucleus, suggesting that mutation hotspots promote increased mutation targeting of spatially co-localized loci (termed hypermutation by proxy). Only around 1% of somatic small variants were in protein coding sequences, but in about half of the driver genes, a contribution of B-cell specific mutational processes to their mutations was found. The B-cell-specific mutational processes contribute to both lymphoma initiation and intratumoral heterogeneity. Overall, we demonstrate that mutational processes involved in the development of gcBCL are more complex than previously appreciated, and that B cell-specific mutational processes contribute via diverse mechanisms to lymphomagenesis.

Fig. 1. Mutation density and replication timing.

A Rainfall plots of three samples including one FL (uppermost track) and two DLBCLs (second and third tracks from the top). For every track, the x-axis displays the genomic coordinate and the y-axis the log-scaled intermutation distance. Clusters of hypermutation (kataegis clusters) can be identified as “rainfalls” reaching very low intermutation distance. The IG loci are highlighted by red vertical lines and red labels, some hallmark genes involved in lymphomagenesis are highlighted by black vertical lines and black labels. B–D Correlation with replication timing. Replication timing is indicated as RepliSeq score of the respective genomic region as determined in [24] (see “Methods” for details). B Scatterplot of replication timing vs. mutation density, showing an inverse relationship between these two quantities. Outliers in this plot, i.e., exceptions from the inverse relationship, are typical targets of SHM in gcBCL. C Boxplot and violin plot of replication time vs. cluster category, demonstrating that kataegis is significantly enriched in early replicating regions (p < 10⁻¹⁶) and psichales in late replicating regions (p < 10⁻¹⁶) of the genome. D Rewiring of replication timing: kataegis regions are located in regions of the genome which are earlier replicating in lymphoblastoid cell lines (y-axis) than in other cell lines (HeLa-S3 (cervical adenocarcinoma), HUVEC (umbilical vein endothelial cells), K562 (chronic myelogenous leukemia in blast crisis), NHEK (epidermal keratinocytes), MCF-7 (mammary gland, adenocarcinoma), IMR-90 (fetal lung fibroblasts), and HepG2 (hepatocellular carcinoma) (x-axis). The light blue color in the background indicates the 95% quantile, the dark blue one the 68% quantile (respective fractions of all SNVs are situated on the colored areas). Regions with a difference in RepliSeq score > 3 are annotated by the closest gene.

Fig. 2. Analysis of mutation density dissects aberrant targeting of SHM and CSR.

A Patterns of nucleotide exchanges in their triplet contexts as extracted cohort wide in the switch regions (upper track) and the regions containing V, D and J genes (middle track). These patterns are not mutational signatures, instead they correspond to visualizations of mutational catalogs. Scales on the y-axes in the different tracks are not fixed, instead a horizontal line is inserted at 5% for rough orientation and comparison. B Clustering of the kataegis clusters according to their contributions from CSR-like and SHM-like mutational processes with contributions of SHM-like and CSR-like as axes. Assessment of the contributions of these two mechanisms to all kataegis clusters was performed by non-negative least squares and subsequent unsupervised k-means clustering (k = 3). Kataegis clusters dominated by a CSR-like pattern are colored in orange, clusters dominated by a SHM-like pattern are colored in green and clusters dominated by neither pattern (other) are colored in purple. C kataegis clusters and kataegis regions displayed as oncoprint. The x-axis encodes samples, the y-axis the kataegis regions, which are ordered by recurrency of affection (≥3%, note that for a better overview, the well established kataegis regions in the IG, BCL2 and BCL6 loci are excluded from the inferred oncoprint-like ordering of the samples and only shown for completeness in the lowest five rows). The oncoprint carries four layers of annotation (normalized horizontal stacked barplots): (i) the fractions of the different kataegis cluster categories (SHM-like = green, CSR-like = orange and other = purple); (ii) the mean distance to the closest TSS in bp; (iii) the fraction of variants overlapping exons (black); and (iv) the fractions of chromatin states from GC B cells annotated to the variants in the respective kataegis regions.

Fig. 3. Hypermutation by proxy (HbP).

A Genome-wide circos diagram showing the positions of all kataegis clusters and their co-occurrence by red arcs. The transparency of these arcs encodes the recurrency of co-occurrence. Arcs are directed from the subject (i.e., primary target) to the object (i.e., secondary target) of the HbP relationship. B–D Detailed illustration of the HbP relationship between S1PR2 and DNMT1. B Co-occurrence: black squares indicate in which samples kataegis clusters are present. Annotation data shows which subgroup the samples belong to, which cell of origin they have and whether a SV is present (DEL_subjectObject: deletion involving both subject (in this case S1PR2) and object (in this case DNMT1); DEL_subject: deletion involving only the subject; TRA_BPsubject: translocation with breakpoint in the subject; DUP_BPsubject: duplication with breakpoint in the subject). C Co-expression in the different subgroups (and normal B cells, other B cells standing for naïve B cells) and D tandem RNA chimeras as detected from RNA-seq: tracks displaying from top to bottom: i) known transcripts of S1PR2 and DNMT1; and Sashimi plots for transcriptomic data of ii) normal GC B cells; iii) lymphoma samples with only S1PR2, i.e., the subject, affected by kataegis; iv) lymphoma samples with both S1PR2 and DNMT1, i.e., subject and object, affected by kataegis; v) lymphoma samples with only DNMT1, i.e., the object, affected by kataegis; vi) lymphoma samples with a deletion affecting either kataegis regions; vii) lymphoma samples with a duplication affecting either kataegis region; and viii) lymphoma samples affected by no event at all in this genomic region. Vertical shading highlights the genomic positions of the two kataegis regions. Numbers on arcs in the sashimi plots display the mean number of splice events (spliced reads) found in the corresponding group.

Fig. 4. New mutational signatures are partially linked to B-cell-specifc mutagenic effects and exhibit characteristic enrichment and depletion patterns.

A Absolute exposures of the samples to the mutational signatures extracted from the combined supervised and unsupervised analyses of mutational signatures. Heights of the stacked bar plots correspond to the number of SNVs explained by the respective mutational signatures. Samples are ordered by subgroup and then decreasing mutational load. For explanation of the identified mutational signatures please refer to the main text. (B, insert) 96-dimensional vectors of nucleotide exchange patterns in the triplet context for the mutational signatures AC9, L1, L2 (all of which were related to AID activity) and L3. Scales on the y-axes in the different tracks are not fixed, instead a horizontal line is inserted at 5% for rough orientation and comparison. C–H Enrichment and depletion patterns of mutational signatures by stratified analyses along different stratification axes, where different colors represent the different strata. C Stratification by genomic regions in which the SNVs were located (“none” = gray – outside of the IG loci, “IG_VDJ” = red – in VDJ genes or intergenic regions between these, “IG_const_switch” = blue – in the switch regions defined in this work, “IG_const_noSwitch” = light blue – in the constant domain of IGH, but outside of the switch regions). Signature L1 is enriched in the switch regions, L2 in the VDJ regions. AC9 is enriched in the constant, non-switch regions (p_KW = 2.2 × 10⁻¹¹, p_Nem = 4.5 × 10⁻⁶). D Stratification by annotated GC B cell-specific chromatin state. L1 was enriched in promoters (p_KW = 1.2 × 10⁻¹⁵, p_Nem = 5.9 × 10⁻¹⁴), while L2 was enriched in transcribed regions (p_KW = 7.7 × 10⁻³⁰, p_Nem = 4.7 × 10⁻¹⁴) and enhancers (p_Nem = 1.4 × 10⁻⁸) as compared to heterochromatic regions. E Stratification by replication timing, illustrating a rewiring of this measure: L1 showed a strong (p_KW = 2.7 × 10⁻¹⁹, p_Nem = 5.3 × 10⁻¹⁴, fold change FC = 1.606) and L2 a moderate (p_KW = 5.7 × 10⁻¹¹, p_Nem = 5.6 × 10⁻⁶, FC = 1.347) enrichment in early replicating regions, as opposed to AC9 which is enriched in late replicating regions (p_KW = 6.3 × 10⁻⁵¹, p_Nem < 2 × 10⁻¹⁶). RS: RepliSeq score. F Enrichment and depletion patterns by subgroup of gcBCL: FLs had higher contributions of L1 (p_KW = 4.74 × 10⁻³, p_Nem = 1.2 × 10⁻³), L2 (p_KW = 6.51 × 10⁻⁴, p_Nem = 1.7 × 10⁻⁴) and AC1 (p_KW = 1.21 × 10⁻⁵, p_Nem = 1.3 × 10⁻⁶) but lower contributions of AC17 (p_KW = 2.02 × 10⁻³, p_Nem = 1.1 × 10⁻³), AC10 (p_KW = 6.51 × 10⁻³, p_Nem = 2.1 × 10⁻⁴), AC6 (p = 2.05 × 10⁻²) and AC2 (p_KW = 8.93 × 10⁻³, p_Nem = 1.3 × 10⁻²) as compared to DLBCLs. G Stratification by consensus clustering of the whole gcBCL cohort. While L3 (p_KW = 9.78 × 10⁻³, enriched in the SOCS1-like, B2M-like and TP53-like consensus clusters), AC1 (p_KW = 9.78 × 10⁻³, enriched in the CSMD1-like and BCL2-like consensus clusters) and AC2 (p_KW = 4.88 × 10⁻², depleted in the BCL2-like cluster) were significantly enriched or depleted between the consensus clusters, L1 (high in the PIM1-like, BCL2-like and MYD88-like consensus clusters) and L2 (high in the B2M-like, BCL2-like and CSMD1-like consensus clusters) only showed trends. H Stratification by consensus clustering of only the DLBCL subgroup. After correcting for multiple testing, no significant effect was observed, with trends for L3 (high in TP53-like) and AC1 (high in BCL2-like). Error bars in C–H display standard error of the mean (SEM).

Fig. 5. B cell-specific mutagenesis alone is not sufficient to drive lymphomagenesis.

Oncoprint of coding (upper part of the figure) and noncoding (lower fifth of the figure dominated by blue color) mutations. The x-coordinate encodes samples which are pre-sorted by subgroups. The y-coordinate encodes different genes or non-coding genes. Different mutation types are encoded by the fill color of the fields in the oncoprint, where different types of mutation can coexist in one sample. Four layers of annotation on the right side of the oncoprint display (i) whether a gene is identified as a driver and (ii) how strongly mutations in AID-specific motifs are enriched, (iii) the best matching signature, and (iv) replication timing.