Multi-omics profiling of longitudinal samples reveals early genomic changes in follicular lymphoma

Abstract

Follicular lymphoma (FL) is the most common indolent type of B-cell non-Hodgkin lymphoma. Advances in treatment have improved overall survival, but early relapse or transformation to aggressive disease is associated with inferior outcome. To identify early genetic events and track tumor clonal evolution, we performed multi-omics analysis of 94 longitudinal biopsies from 44 FL patients; 22 with transformation (tFL) and 22 with relapse without transformation (nFL). Deep whole-exome sequencing confirmed recurrent mutations in genes encoding epigenetic regulators (CREBBP, KMT2D, EZH2, EP300), with similar mutational landscape in nFL and tFL patients. Calculation of genomic distances between longitudinal samples revealed complex evolutionary patterns in both subgroups. CREBBP and KMT2D mutations were identified as genetic events that occur early in the disease course, and cases with CREBBP KAT domain mutations had low risk of transformation. Gains in chromosomes 12 and 18 (TCF4), and loss in 6q were identified as early and stable copy number alterations. Identification of such early and stable genetic events may provide opportunities for early disease detection and disease monitoring. Integrative analysis revealed that tumors with EZH2 mutations exhibited reduced gene expression of numerous histone genes, including histone linker genes. This might contribute to the epigenetic dysregulation in FL.

Fig. 1. Landscape of genetic alterations in nFL and tFL.

A Schematic overview of the FL multi-omics data cohort. The FL cohort was enriched for patients who experienced transformation (tFL group) as compared to FL patients who relapsed without transformation (nFL group). The longitudinal samples were subjected to multi-omics analysis by WES, SNP array, and RNA-seq. B Recurrent and significantly mutated genes in nFL and tFL patients as identified from WES (n = 94 samples) by at least one of the three applied tools (MutSig2CV, IntOgen pipeline and 2020Plus). Only cancer driver genes appearing at average VAF > 0.15 (adjusted for tumor content) are shown. The genes are ordered based on their mutational frequency. Mutation types are color-coded as indicated; for genes with multiple mutations, the most severe mutation with consequence for change in the protein structure is shown. The FL classification was performed as described in Dreval et al. [28]. C Forest plots of cancer driver gene mutations with odds ratio for association with nFL group vs. tFL group; n = 44, statistical testing by Fisher exact test. Left: the unified mutations per patient were used for the test. Right: the mutations in the pretreatment biopsies were used for the test. Mutant EZH2 was significantly enriched in the tFL group. D-E Kaplan–Meier curves showing D overall survival and E time to transformation between cases with CREBBP KAT domain mutation vs. the other cases in the WES discovery cohort.

Fig. 2. Clonal shifts are common during FL evolution as demonstrated by large genomic distances between biopsies in individual FL patients.

Tumor genome evolution was inferred with PyClone from ultra-deep WES of longitudinal samples, and clonal phylogenies were constructed for all tumors from the same patient, based on the clustering patterns and cluster means using ClonEvol. Genomic distances were calculated on the variants with cellular fraction above 0.25 as inferred by PyClone. A The nFL patient P21 and tFL patient P32 are examples of divergent evolution: The two biopsies of P21 are represented by disjoint sets of clones. The subclonal composition of all three biopsies from tFL patient P32 suggests divergent evolution from a CPC as the most prevalent clone of each biopsy was absent from the other biopsies. B The nFL patient P14 and the tFL patient P29 are examples of linear evolution. For P14, the relapsed tumor evolved directly from the major clone in the pretreatment biopsy (clone 1). For P29, the two biopsies shared the dominant clone (clone 1). C Genomic distances between any two biopsies of each patient were calculated on the variants with cellular fraction above 0.25 as inferred by PyClone. Shown is genomic distance vs. time interval of the compared serial biopsies and D Genomic distance ordered on the x-axis based on the largest genomic distance between any two biopsies of a patient. E Comparison of the genomic distances in nFL cases (left) or tFL cases (right) classified as divergent or linear evolution. F Comparison of the genomic distances when comparing nFL vs tFL patients (left), or comparing two FL biopsies vs. FL and DLBCL biopsies within individual patients (middle), or when comparing biopsies from non-POD24 patients vs. biopsies from POD24 patients (right). Statistical difference by Student t-test, p < 0.05.

Fig. 3. CREBBP and KMT2D are early dominant clonal mutations in FL.

Tumor genome evolution was inferred with PyClone and CloneEvol as in Fig. 2. A Illustration of the mapping of Dominant Clonal Mutations (DCMs) to divergent phylogenetic trees from a common progenitor cell (CPC). Mutations found in all tumor clones, such as “A”, are assumed to have arisen in CPC, and were defined as early DCMs. Other mutations found in one or more (but not all) dominant clones were defined as late DCMs, such as “B” and “C”. B Genes harboring dominant clone mutations, DCM genes. DCM genes mutated in at least four patients out of 32 are shown.

Fig. 4. Identification of early and stable copy number changes in FL.

A Bubble plot identifying all genomic regions, by genomic location, in which an amplification was identified at the first available biopsy and the copy number state remained within ±0.3 logR or alterations above 0.3 logR increased in magnitude across all available biopsies of a patient. B Bubble plot identifying all genomic regions, plotted by genomic location, in which a deletion was identified at the first available biopsy and the copy state stayed within ±0.3 logR or alterations below −0.3 logR increased in magnitude across all available biopsies of a patient. The y-axis represents the number of cases in which that genomic region was identified as an early and stable copy number region. For A–B, the bubble plot size represents the size of the genomic region identified and color indicates the ratio of tFL cases over nFL cases identified as containing the region as an early stable copy state. SNP6.0 data was available for 34 patients with serial biopsies. C 2D plot showing the Pearson correlations for cis-genes in the discovery and the Nordic validation cohort. In addition, significantly differentially expressed cis-genes with a Pearson correlation above 0.4 in both cohorts were highlighted in orange. D–F Correlation of omics data and gene dependency data for TCF4, using DepMap portal [68]. D TCF4 RNA expression and protein expression. E TCF4 RNA expression with RNA interference data, showing dependency of cell lines on gene presence. Negative values below −0.5 correlate with depletion of TCF4 with decreasing cell proliferation rates. F TCF4 RNA expression with CRISPR dependency scores (Chronos). Chronos scores below −0.5 correlate with depletion of TCF4 with decreasing cell proliferation rates. For E-F, the linear regression was correlated for the 3 highlighted lineage types, the R squared and p-value for the linear fit is calculated.

Fig. 5. Histone genes are downregulated in FL.

Association between EZH2 mutation and downregulation of histone genes. B-cell transcriptomes were inferred by CIBERSORTx from RNA-seq data (this study, n = 64) and from the Nordic cohort with gene expression data described in Steen et al. [27] (n = 82). Volcano plot displaying significant differentially expressed genes (DEG) in cases with EZH2 mutations vs. wild-type in A left plot: the RNA-seq cohort (this study), middle plot: the Nordic cohort [27]. In the right plot, analysis of scRNA-seq data from Han et al. [45] is showing DEG in malignant B cells from EZH2 mutant cases vs. wild-type cases. Significant DEG were defined as −0.35 < FC(log₂)>0.35 and p < 0.05. B Unsupervised hierarchical clustering analysis of the expression level of histone genes in the RNA-seq cohort of FL patients (this study) and in FACS-sorted B-cell populations from healthy donor tonsils. Three clusters with generally higher (HL-H), intermediate (HL-M) and lower (HL-L) histone gene expression patterns were identified.