Prognostic gene expression and microRNA profiling signatures and genetic alterations in primary testicular diffuse large B-cell lymphoma

Abstract

Primary testicular (PT) diffuse large B-cell lymphoma (DLBCL) is a rare and aggressive lymphoma with distinct clinical and molecular characteristics. To identify prognostic biomarkers in PT-DLBCL, in this study we analyzed DNA and RNA samples of PT-DLBCL tumors from 206 patients using next-generation sequencing platforms and assays. Genetic alteration analysis found that multiple chromosomal copy number variations (CNVs), TP53 transcript mutations with high variant allele frequency, and MCD subtype had significantly adverse prognostic effects, whereas elevated microsatellite instability had a significantly favorable prognostic effect in PT-DLBCL. Targeted RNA-seq analysis identified a PTL gene expression signature by comparing PT-DLBCL with systemic DLBCL and revealed the heterogeneity within PT-DLBCL by unsupervised clustering, which classified PT-DLBCLs into a testicular lymphoma tumor (TLT) subtype and a microenvironment (ME) subtype. The TLT subtype featured upregulation of genes functioning in DNA damage response, DNA repair, chromatin remodeling, the cell cycle, and the nucleus, and was associated with significantly poorer patient survival and higher frequencies of MYD88 mutations, multiple CNVs, MCD subtype, bulk tumors, and elderly patients in the PT-DLBCL cohort. In contrast, the ME subtype distinctively featured upregulation of various signaling pathway genes involving the tumor microenvironment and downregulation of BTK and B-cell receptor signaling genes, and was associated with significantly better clinical outcome than the TLT subtype of PT-DLBCL independently of CNVs, MCD and MYD88 mutation and than systemic DLBCL. Moreover, genomic microRNA profiling analysis identified a PTL microRNA signature significantly differentially expressed between PT-DLBCL and systemic DLBCL patients and within the PT-DLBCL cohort, and PT-DLBCL patients with higher expression of 16 PTL microRNAs (14 are testicular tissue-specific) had significantly better survival. In summary, this study revealed the molecular heterogeneity in genetic abnormalities and expression profiles of coding genes and microRNAs within the PT-DLBCL entity, and identified significant prognostic biomarkers and PTL signatures.

Fig. 1. Genetic alteration analysis and identified biomarkers in PT-DLBCL.

A Mutation overview and the most frequently mutated genes by targeted DNA-sequencing in 192 patients with PT-DLBCL. The VAF cutoff was 5%. B Oncoplot showing the case distribution of LymphGen subtypes and some genetic hallmarks in 192 patients with PT-DLBCL. The most frequent MCD genetic subtype was associated with significantly poorer overall survival (OS) in the PT-DLBCL cohort. Note: VAF cutoff was >2% for LymphGen subtyping; compared with results in panel A, 63 mutations (in 27 genes) with a VAF 2.05-4.92% were not excluded for genetic subtyping. C Remarkably, presence of TP53 RNA variants with a VAF ≥ 40% by RNA-seq analysis was associated with significantly poorer OS and progression-free survival (PFS). D Patients with chromosome 12 gain, chromosome 19 copy number alterations (CNVs), or ≥3 chromosomes with CNVs had significantly poorer OS rates. E Scatter plot: PT-DLBCL patients with a high VAF of TP53 RNA variants by RNA-seq analysis had significantly higher mean and median (by 2-tailed t test and Mann-Whitney (U) test, respectively) CNV numbers in DNA-sequencing analysis. Bars show mean values, lines show ±SD. The first violin plot: PT-DLBCL patients with TP53 variants by RNA-seq analysis (VAF all ≥10%) had significantly higher mean and median numbers of CNVs in DNA-sequencing analysis. The second violin plot: PT-DLBCL patients with ≥3 chromosomal CNVs were associated with lower proportions of CD8⁺ T cells and M0 macrophages and higher proportions of naïve B cells, resting memory CD4 T cells, and activated NK cells, quantified through RNA-seq data deconvolution by CIBERSORT (refs. [40, 41]). Lines in violin plots represent medians and quantiles in Mann-Whitney (U) test. All P values are 2-tailed. F A combined group of patients with BCL2 gene amplification, >4 polyploid, or BCL2 gene rearrangement by FISH analysis had significantly poorer OS than other patients in the PT-DLBCL cohort. G High percentages of MSI sites among evaluated sites (>100 in most cases) in DNA-sequencing analysis were associated with significantly better OS of patients with PT-DLBCL.

Fig. 2. Targeted gene expression profiling and prognostic analysis in PT-DLBCL.

A Kaplan-Meier survival curves from left to right: The sequenced PT-DLBCL and systemic DLBCL cohorts did not show a significant difference in overall survival (OS). The ME subtype of PT-DLBCL had significantly better OS than the TLT subtype of PT-DLBCL and systemic DLBCL cases. When both age and GEP subtypes were integrated into the survival analysis, the favorable prognostic effect of the ME subtype was significant only in patients ≤60 years old. When both MCD (DNA-sequencing subtype) and TLT/ME subtypes (RNA-sequencing subtypes) were integrated into the survival analysis, the TLT subtype showed a significantly adverse prognostic impact in MCD patients; the MCD subtype was associated with a statistically non-significant trend of poorer OS in the TLT subtype. B Left heatmap: expression of significantly differentially expressed genes (DEGs) between PT-DLBCL and systemic DLBCL patients. Right heatmap: median-centered unsupervised clustering based on expression of the 150-gene PTL signature in the PT-DLBCL cohort formed two clusters of patients (each column represents one patient), as well as two clusters of genes, 45 genes (TLT subsignature) and 105 genes (ME subsignature) of the PTL signature (each row represents expression of one gene), respectively. The vertical red/green color bars for two gene clusters on the left side of this heatmap and the lighter color bars on the right side of the heatmap in the left panel indicate that the gene clusters were similar but a few genes were different. The case plot aligned underneath the heatmap shows the status of MYD88 mutation detected by DNA and RNA sequencing, respectively, TP53 mutation detected by RNA sequencing with VAF ≥ 10% and VAF ≥ 40%, respectively, MSI-high status, ≥3 chromosomal CNVs by DNA-sequencing analysis, and LymphGen genetic subtypes for patients in the heatmap. Sequencing-failed cases for MYD88/TP53 mutations, MSI-high, and ≥3 CNVs status were in grey color, and were in white color for LymphGen genetic subtypes. Significance of frequency differences between TLT and ME patient clusters: *P < 0.05; ***P < 0.001; ****P < 0.0001. C Enrichment dot plots to visualize part of the significantly enriched GO and KW terms and KEGG pathways (FDR < 0.05) by DAVID functional annotation for 45 genes (TLT signature), 105 genes (ME signature), 384 genes (Mut-MYD88 signature in PT-DLBCL), and 121 genes (downregulated in cases with ≥3 CNVs. Supplementary Tables 1, 3, and 4). The color scale shows the enrichment P values, the horizontal axis shows fold enrichment or gene ratio, and dot sizes represent numbers of genes hits for the significant GO and KW terms. GO Gene Ontology, KW UniProtKB Keywords, KEGG Kyoto Encyclopedia of Genes and Genomes, BP biological process, CC cell component, MF molecular function, EASE Expression Analysis Systematic Explorer (A desktop version of DAVID). D Left heatmap: expression of 447 DEGs between patients with and without MYD88 RNA mutations. Right heatmap: expression of 134 DEGs between patients with and without ≥3 CNVs by DNA-sequencing analysis.

Fig. 3. GEP CIBERSORT, clustering, and prognostic analysis.

A Violin plots showing normalized BTK and IL15 levels in the heatmaps in Fig. 2B. The ME subtype had significantly lower BTK transcription than the TLT subtype of PT-DLBCL and systemic DLBCL (sys-DLBCL). The TLT subtype had significantly lower IL15 transcription than the ME subtype of PT-DLBCL and systemic DLBCL. Lines represent medians and quantiles, P values are by 2-tailed Mann-Whitney (U) test. B Scatter plot showing significant differences in CIBERSORT leukocyte deconvolution results between the TLT and ME clusters. Each dot represents one case; bars and lines represent mean proportions ±SD. Significance: *P < 0.05; **P < 0.01; ***P < 0.001, ****P < 0.0001. C Left heatmap: expression of 57 genes upregulated and 109 genes downregulated in PT-DLBCL vs. systemic DLBCL with MYD88 mutations (Supplementary Table 7). The Venn diagram shows the overlaps of these DEGs with the PTL signature. Right heatmap: median-centered unsupervised clustering based on expression of 57 genes formed a distinct cluster with higher expression of 21 genes included in the ME-PTL subsignature. This cluster was associated with significantly better patient survival in the PT-DLBCL cohort. Similar unsupervised clustering was performed in the Mut-MYD88 PT-DLBCL subcohort (Supplementary Fig. 6A), and the resultant cluster with high expression of ME subsignature was associated with significantly better patient survival as shown by the Kaplan-Meier survival curves. D Consensus clustering in MCD-like PT-DLBCL (subtyped by only hallmark mutations and BCL2/BCL6/MYC fusions but not by CNVs, Supplementary Fig. 6B). Left: consensus matrix heatmap based on the optimal clusters value (k = 5) by cophenetic correlation coefficients for hierarchically clustered matrices. Middle: enriched pathways by Gene Set Variation Analysis in the five clusters/groups of patients. Significant enrichment by multiple correction was marked by an asterisk. Right: survival curves for patients in the five clusters/groups.

Fig. 4. MicroRNA profiling analysis.

A Heatmap: significantly differentially expressed 99 microRNAs between profiled systemic DLBCL and PT-DLBCL cases. Boxplots: top 10 upregulated and top 10 downregulated microRNAs in PT-DLBCL vs. systemic DLBCL. B Heatmap showing top 40 significantly differentially expressed microRNAs (including 30 up and 10 down microRNAs) between systemic DLBCL and PT-DLBCL. C Heatmap: median-centered unsupervised clustering of sequenced PT-DLBCL cases based on expression of 35 microRNAs (including 16 microRNAs up- and 19 down-regulated in PT-DLBCL). Kaplan-Meier survival curves: patients with high expression of PT-DLBCL-associated 16-microRNA signature had significantly better survival. D Heatmap for median-centered unsupervised clustering of sequenced PT-DLBCL and systemic DLBCL cases using 35 microRNAs significantly differentially expressed between them. E Read counts of 16 microRNAs in samples of testis, blood, and lymph nodes in the miTED database.