The Mutation of BTG2 Gene Predicts a Poor Outcome in Primary Testicular Diffuse Large B-Cell Lymphoma.

Introduction: Primary testicular diffuse large B-cell lymphoma (PT-DLBCL) is a rare and aggressive form of mature B-cell lymphoma commonly found in elder males, but its genetic features are poorly understood. In this study, we had performed target-sequencing of 360 lymphoma-related genes on 76 PT-DLBCL patients with a median age of 65 (33-89). Our data provide a comprehensive understanding of the landscape of mutations in a small subset of PT-DLBCL.
Methods: A total of 76 PT-DLBCL patients were sequenced, and their clinical data and follow-up data were collected. The relationship between mutated genes, clinical data and prognosis and survival of PT-DLBCL patients was retrospectively analyzed by statistical software.
Results: We observed a median of 15 protein-altering variants per patient in our data and was identified recurrent oncogenic mutations of 360 lymphoma-related genes involved in PT-DLBCL, including PIM1 (74%), MYD88 (50%), KMT2D (38%), KMT2C (34%), BTG2 (34%), TBL1XR1 (34%) and ETV6 (24%). Compared with classic DLBCL, PT-DLBCL showed an increased mutation frequency of PIM1, MYD88, BTG2, while NOTCH1 appeared exclusive mutated with PIM1, MSH3 and ETV6. Cox risk model regression analysis showed that age ≥60 years, IPI 3-5 points, BTG2 gene mutation and extranodal organ invasion suggested poor prognosis. Finally, we constructed an OS predict model of PT-DLBCL patients using above factors with a high accuracy.
Conclusion: In conclusion, our results revealed genomic characterization of PT-DLBCL, and the mutation of BTG2 was an independent factor predicting a poor prognosis.

Introduction

Primary testicular diffuse large B-cell lymphoma (PT-DLBCL) is a rare and aggressive form of mature B-cell lymphoma.1–3 PT-DLBCL was the most common type of testicular tumor in men aged over 60 and characterized by painless uni- or bilateral testicular masses with infrequent constitutional symptoms.4–6 PT-DLBCL shows significant extranodal tropism, and features a high risk of recurrence in the central nervous system (CNS) and contralateral testis, directly leading to poor prognosis.7,8 Although the addition of radiotherapy, whole-course chemotherapy, CNS-directed prophylaxis and rituximab have greatly improved the prognosis of DLBCL patients,9–11 the biomarkers of prognosis remains poor for PT-DLBCL.

Gene expression profiling had demonstrated that PT-DLBCL showed phenotypic similarities with nodal DLBCL of activated B-cell-like (ABC) origin.12–14 Research based on small samples found that mutations of TP53 and C-MYC as well as BCL2 rearrangements were rare in PT-DLBCL, whereas BCL6 was commonly rearranged and CXCR4 might be the prognostic marker for PT-DLBCL.15,16 Chapuy et al have found that although patients with PT-DLBCL had few TP53 mutations compared with systemic DLBCL ones, but the P53 pathway had always been interfered via disorders of upstream signals, usually double-allelic CDKN2A deletion.17 Moreover, patients with central nervous system lymphoma (PCNSL) and primary testicular lymphoma (PTL) were observed with high-frequent mutation of CD79B, MYD88L265P, ETV6, PIM1 and TBL1XR1 mutations, as well as poor prognosis.18,19 However, due to the rarity of PT-DLBCL (approximately 1–2% of cases and rarity of specimen),20 the genetic features of PT-DLBCL had not been understood well, and the effects of gene mutation on the prognosis were still unclear.

Previous studies have shown that type B symptoms, advanced neuronal cell necrosis (III/IV), and extra node involvement were poor prognostic markers for PT-DLBCL.21–26 Retrospective studies had shown that the prognostic roles of IPI may be surrogate indicators of high tumor load and diffuse disease. For most of PT-DLBCL patients in early stage, IPI score of them (always equated or less than 2) had limited prognostic value. In addition, mutations of MYD88 and CD79B were frequently observed in PT-DLBCL, but the prognostic effects of them were not rarely reported.27

Herein, in the present study, we performed targeted sequencing of 360 lymphoma-related genes on 76 patients with PT-DLBCL and analyzed their mutation profiles with clinical indicators, immunohistochemical, and prognostic outcomes. Furthermore, we constructed of a prognostic model (PI, riskScore) based on clinical indicators and gene mutation. Our study could be benefit for the understanding of mutation profiles of lymphoma-related genes in PT-DLBCL and provide valuable therapeutic targets.

Materials and Methods

Data Source

A total of 76 PT-DLBCL patients admitted to the Affiliated Hospital of Nantong University and Tumor Hospital Affiliated to Nantong University from January 2007 to December 2020 were selected retrospectively as study subjects. The study was approved by the Ethics Committee of the Affiliated Hospital of Nantong University and Tumor Hospital Affiliated to Nantong University in accordance with the Helsinki Declaration (Reference number: 2020-L105; 2021-081). All the patients selected for our study were fully informed about our protocols and signed an informed consent to participate in this study.

All of the patients enrolled in this study were diagnosed by histopathology and immunohistochemistry in the pathology department of the Affiliated Hospital of Nantong University and Tumor Hospital Affiliated to Nantong University. The histopathological characteristics of each sample had been represented in the . The median age of patients was 65 (33–89), and most of patients underwent orchiectomy, followed by immunochemotherapy (R-CHOP), contralateral testicular irradiation (25–30Gy) and prevention of central nervous system. The median follow-up time of 39 (1–124) months.

Library Construction

DNA samples were extracted from pathological biopsies of testicular tissues. DNA samples were fragmented via Bioruptor (Diagenode, Bioruptor UCD-200) following manufacturer’s protocol. Libraries were constructed using the KAPA Hyper DNA Library Prep Kit (KAPA Biosystems, KK8504). At last dual-indexed sequencing libraries were PCR amplified with KAPA HiFi Hot start-ready Mix (KAPA, KK2602) for 4–6 cycles, then cleaned up by purification Beads (Corning, AxyPrep Fragment Select-I kit, 14223162). Library concentration and quality were determined by the Qubit 3.0 system (Invitrogen) and Bioanalyzer 2100 (Agilent, Agilent HS DNA Reagent, 5067-4627).

Hybrid Selection and Ultra-Deep Next Generation Sequencing

The 5′-biotinylated probe solution was used as the capture probes. The probes for targeted sequencing cover exons and selected introns of 360 lymphoma-related genes (). Briefly, 1 μg of each fragment sequencing library was mixed with 5 μg of salmon sperm DNA, 5 μg of human Cot-1 DNA, and 1 unit adaptor-specific blocker DNA in hybridization buffer, heated for 10 minutes at 95°C, and held for 5 minutes at 65°C in the thermocycler. The capture probes were added to the mixture within 5 minutes, and the solution hybridization was performed for 16–18 hours at 65°C. After hybridization, the captured targets were selected by pulling down the biotinylated probe/target hybrids using streptavidin-coated magnetic beads, and off-target library was removed using wash buffer. The PCR master mix was directly added to amplify (6–8 cycles) the captured library from the washed beads. After that, the samples were purified by AMPure XP beads, quantified by qPCR (Kapa) and sized on bioanalyzer 2100 (Agilent, Agilent HS DNA Reagent, 5067-4627). Libraries were normalized to 2.5 nM and pooled. Finally, the library was sequenced as paired 150-bp reads on Illumina HiSeq 4000 according to manufacturer’s instrument. The sequencing depth of genes varied from 20 to 6821 between patients, and the detailed values of sequencing depth were deposited in the .

Sequence Alignment and Processing

Base calling was performed on bcl2fastq v2.16.0.10 (Illumina, Inc.) to generate sequence reads in FASTQ format (Illumina 1.8+ encoding). Quality control (QC) was performed with Trimmomatic. High quality reads were mapped to the human genome (hg19, GRCh37 Genome Reference Consortium Human Reference 37) using the BWA aligner 0.7.12 with BWA-MEM algorithm and default parameters to create SAM files. Picard 1.119 was used to convert SAM files to compressed BAM files which were then sorted according to chromosome coordinates. The Genome Analysis Toolkit (GATK, version 3.4-0) was used to locally realign the BAMs files at intervals with indel mismatches and recalibrate base quality scores of reads in BAM files.

SNVs/Indels/CNVs Detection

Single nucleotide variants (SNVs) and short insertions/deletions (indels) were identified by VarScan2 2.3.9 with minimum variant allele frequency threshold set at 0.01, and p-value threshold for calling variants set at 0.05 to generate Variant Call Format (VCF) files. All SNVs/indels were annotated with ANNOVAR, and each SNV/indel was manually checked on the Integrative Genomics Viewer (IGV). Copy number variations (CNVs) were detected using in-house-developed software. SNV re-filtering criteria:1. ≥5 reads and variant allele fraction (VAF) of ≥1%, 2. 1000 Genomes Project or ExAC database >1% filtered out. INDEL re-filtering criteria: 1. ≥4 reads and variant allele fraction (VAF) of ≥1%; 2. 1000 Genomes Project or ExAC database >1% filtered out. CNV re-filtration criteria: 1. amplification factor ≥ 2; 2. single copy number deletion ≤ 0.6; 3. double copy deletion ≤ 0.2.

SVs Detections

We used Delly fusion calling tool to identify the number of chimeric reads (sequencing paired ends mapped to different genes) and split reads (spanning a fusion breakpoint) from the targeted DNA-seq data.

Statistical Analysis

Statistical methods were all analyzed by SPSS 26.0 and GraphPad 8.0. P values were all two-sided and statistical significance was set at P < 0.05 if not mentioned. All confidence intervals were stated at the 95% confidence level. Mann–Whitney U-test were used for categorical variables to compare differences between groups, and nonparametric rank sum test was used for continuous variables. The overall survival (OS, Overall survival) was analyzed and plotted by Kaplan-Meier method, and the Log rank test was used to compare the survival rates among the groups. Cox proportional hazards model was used to estimate the HR and its associated confidence interval (CI).

Results

Patient Characteristics

A total of 76 PT-DLBCL patients with a median age of 65 (range 33–86) years were included in this study. According to the Ann Arbor staging system, 43 patients and 33 patients were diagnosed with stage I/II and stage III/IV disease, respectively. There were 9 patients (11.8%) who were diagnosed with the germinal center B-cell-like (GCB) and 67 (88.2%) with the non-GCB subtype. According to the IPI score, 40 cases were in the low-risk group and 11 cases were in the low-medium risk group. There were 10 cases in the medium-high risk group and 15 cases in the high-risk group. 24 patients (31.6%) had extranodal lesions. The unilateral testicle was involved in 69 (90.8%) patients, while the bilateral testicular involvement was observed in 7 patients (9.2%). 9 cases were diagnosed with B symptoms (Table 1).

Table 1

Baseline Characteristics of 76 PT -DLBCL Patients

Characteristic	Patients	Percent %
Age
≤60	21	27.6
>60	55	72.4
Ann Arbor stage
I-II	43	56.6
III-IV	33	43.4
Subtypes
GCB	9	11.8
Non-GCB	67	88.2
IPI/score
Low risk	40	52.6
Low-medium risk	11	14.5
Medium-high risk	10	13.2
High risk	15	19.7
Invasion
Yes	24	31.6
No	51	68.4
Location
Unilateral	69	90.8
Bilateral	7	9.2
B symptoms
Yes	9	11.8
No	67	88.2

Mutation Profile of PT-DLBCL

The exon regions of all the genes were amplified. In total, all of the 360 lymphoma-related genes were well sequenced in our samples (Figure 1). We observed that 40 genes were high-frequently mutated in our samples (more than 6 cases, Figure 2). As shown in Figure 3A and B, 75 of 76 (98.68%) samples demonstrated genetic mutations. Missense mutation was the most common variant classification. SNPs were the most common variant type, and C > T ranked as the top SNV class. The results also demonstrated PIM1 (74%) as the gene with the highest mutation frequency, followed by MYD88 (50%), KMT2D (38%), KMT2C (34%), BTG2 (34%), TBL1XR1 (34%), ETV6 (24%), PRDM1 (24%), FAT1 (21%) and DUSP2 (20%). The distribution of gene mutation abundance was displayed by Plotting VAF (Figure 1A) and Rainfall plots showed the hyper mutated genomic region by plotting the mutation distances on each chromosome (Figure 1B). We classified SNPs into transition and transversion, and the ratio of transition exceeded the ratio of transversion in 76 PT-DLBCL patients (Figure 1C).

Figure 1

Visualization of the type and location of gene mutations. (A) Gene mutation abundance distribution of top 10 genes. (B) The hyper mutated genomic region on each chromosome. (C) Differences in frequency distribution of transitions and transversions in patients with PT-DLBCL.

Figure 2

Gene mutation analysis in patients with PT-DLBCL. The number on each bar indicates the number of patients carrying the indicated gene mutation.

Figure 3

(A) Landscape of genetic and expression variation of related genes in PT-DLBCL. (B) The mutation frequency and classification of related genes in PT-DLBCL.

Finally, co-occurrence and exclusivity of the most frequent variants were computed to reveal significant patterns of mutational co-segregation underlying biological cooperativity of certain mutational events (Figure 4). It could be found that co-occurrent variants with statistically significance were TBL1XR1 and BTG2, TBL1XR1 and FOXO1; FAT1 and ARID1A, FAT1 and ETV6, ARID1A and NOTCH1, ETS1 and CIITA. Co-exclusive variants with statistically significance were NOTCH1 and PIM1, MSH3 and ETV6.

Figure 4

Comparison of mutation event correlation of top 20 genes. Co-occurrence or exclusivity of the most recurrent mutational events in the 76 patients were assessed using the Somatic interactions method with a dot plot indicating the odds ratio of co-occurrence (green) or exclusivity (yellow). *P<0.05 and P < 0.1, P<0.05 is statistically significant.

Basic Features of Patients with BTG2 Mutation and BTG2 Wild Type (wt)

From The Cancer Genome Atlas (TCGA) database (), we found the expression distribution of BTG2 gene in DLBCL was high than in normal tissues (Figure 5A). Next, we analyzed the BTG2 mutation of 76 PT-DLBCL samples. BTG2 mutations, found in 26 cases (34%), comprised missense mutations, frameshift deletion and nonsense mutation. The majority of the mutations (84.6%, n=22) were missense variants (Figure 5B). In Figure 5C, protein structure of BTG2 with mutations was illustrated. Mutant residues are highlighted in different color. According to the results of NGS detection, 76 patients were divided into BTG2 mutation group (n=26) and BTG2 wt group (n=50). The basic clinical characteristics were analyzed (Table 2). There was no significant difference between the two groups in age, tumor invasion, disease location, PIM1, MYD88, KMT2C, KMT2D gene mutation rates, but there were significant differences in IPI (P=0.019), pathological type (P=0.021), TBL1XR1 (P=0.009) mutation between the BTG2 mutation group and BTG2 wt group.

Figure 5

Schematic representation of mutations in BTG2. (A) Comparison of BTG2 expression level between DLBCL and normal group. (B) Location of mutations in the BTG2 protein. Frequency of mutations (y-axis) is shown by amino acid (x-axis). (C) Mutant residues of BTG2 are highlighted in different color.

Table 2

Correlation Analysis of Basic Characteristics Between BTG2 Mutation and BTG2 wt Patients

	BTG2 wt	BTG2 Mutation	P value
	N = 50	N = 26
Age
≤60	15	6	0.528
>60	35	20
IPI/score
Low risk	30	10	0.019*
Low-medium risk	9	2
Medium-high risk	3	7
High risk	8	7
Subtypes
GCB	9	0	0.021*
Non-GCB	41	26
Invasion
Yes	14	10	0.359
No	36	16
Location
Unilateral	47	22	0.229
Bilateral	3	4
Mutation
PIM1	33	23	0.065
MYD88	26	12	0.604
KMT2D	19	10	0.969
KMT2C	16	10	0.579
TBL1XR1	12	14	0.009*
DUSP2	19	10	0.969
PRDM1	16	10	0.579
ETV6	14	4	0.409
FAT1	9	7	0.327

Note: *The difference is statistically significant (P < 0.05).

Prognostic Analysis of PT-DLBCL Patients

Univariate analysis of 76 patients who completed the follow-up showed that age > 60 years, Ann Arbor stage III-IV, IPI 3–5 points, extranodal organ invasion, BTG2 mutation, TBL1XR1 mutation and PRDM1 mutation were risk factors for prognosis (P < 0.05) (Table 3). Multivariate Cox proportional hazard regression analysis revealed that BTG2 mutations dominated the death risk of PT-DLBCL patients, with hazard ratios in the range of 95% CI 1.15 to 4.3311 (HR, 2.23, P=0.007). The median OS was 28 months in the BTG2 mutation group and 65 months in the BTG2 wt group (Figure 6A). The prognostic effect of BTG2 mutation on OS was similar (all HR >1.0) across all prespecified subgroups (Figure 6B).

Table 3

Univariate Analysis of OS Influencing Factors in 76 PT-DLBCL Patients

Factors	N	Median OS (Months)	P value
Age
≤60	21	68	0.039*
>60	55	45
Ann Arbor stage
I–II	43	68	0.000*
III–IV	33	29
IPI/score
Low risk (0–1)	40	71	0.000*
Low-medium risk (2)	11	38
Medium-high risk (3)	10	32
High risk (4–5)	15	12
Subtypes
GCB	9	50	0.145
Non-GCB	41	38
Invasion
Yes	24	27	0.001*
No	51	58
Location
Unilateral	69	50	0.366
Bilateral	7	31
PIM1
Mutation	56	50	0.840
Wild type	20	53
MYD88
Mutation	38	50	0.806
Wild type	38	60
BTG2
Mutation	26	28	0.007*
Wild type	50	65
KMT2C
Mutation	26	50	0.361
Wild type	50	68
KMT2D
Mutation	29	50	0.239
Wild type	47	50
TBL1XR1
Mutation	26	28	0.033*
Wild type	50	60
DUSP2
Mutation	15	53	0.862
Wild type	61	50
PRDM1
Mutation	18	48	0.040*
Wild type	58	NA
ETV6
Mutation	18	53	0.862
Wild type	58	50
FAT1
Mutation	16	32	0.494
Wild type	60	52

Note: *The difference is statistically significant (P < 0.05).

Figure 6

Kaplan-Meier Estimates of Overall Survival and Forest Plots by Subgroup. (A) Overall Survival of BTG2 mutation and BTG2 wild type (wt). (B) Overall Survival in Prespecified Subgroups.

Construction of a Prognostic Model Index (PI, riskScore) Based on Clinical Indicators and Gene Mutation

Cox risk model regression analysis showed that age, extranodal organ invasion, IPI score and BTG2 gene mutation were independent factors affecting the prognosis of patients (Table 4). As shown in Figure 7A, the status of BTG2 mutation combined with IPI score can further distinguish the prognosis of PT-DLBCL patients, especially for patients with IPI score<3 (P<0.0001). In addition, a prognostic nomogram was created to quantify the relationship between these risk genes and survival. From this nomogram, we could obtain the total points and estimate the 1-year and 3-year survival rate of each patient (Figure 7B).

Table 4

Cox Risk Model Regression Analysis of OS in 76 PT-DLBCL Patients

Factors	HR	95% CI	P value
Age	2.288	1.058–4.949	0.035*
Invasion	2.113	1.070–4.174	0.031*
IPI/score	1.553	1.225–1.968	0.000*
BTG2	1.883	1.014–3.499	0.045*

Note: *The difference is statistically significant (P < 0.05).

Figure 7

Establishment of a prognostic model. (A) Overall survival of BTG2 mutation combined with IPI score (0, BTG2 wild type and IPI<3; 1, BTG2 mutation and IPI<3; 2, BTG2 wild type and IPI≥3; 3, BTG2 mutation and IPI≥3). The dashed line indicates the median. (B) A nomogram was developed based on the results of multivariate Cox proportional hazards analysis to predict the 1-year and 3-year overall recurrence.

Discussion

Functional mutation and epigenetic alteration of genes are driver events in cancers,28 while the accumulation of necessary somatic mutations is a leading cause of the initiation and development of DLBCL.29 Deep understanding of high-frequent mutated genes could provide clues to better elucidate biological mechanisms of tumorigenesis and identify biomarkers for diagnosis and therapy,30 especially for the rare type of cancers, such as PT-DLBCL. In this study, we comprehensively investigated the mutation profiles of lymphoma-related genes and their impacts on prognosis of PT-DLBCL patients, as well as TCGA DLBCL cohort. Our data demonstrated that the mutation profiles of PT-DLBCL were distinguish from classic DLBCL, and the mutation of BTG2 gene may serve as a valuable biomarker in prognosis of PT-DLBCL.

The functional variation and dysregulation of genes have long been thought associated with malignancies and have been extensively studied in the past decades.31 Our data demonstrated that the high-frequent mutated genes could be relevant to specific pathways, leading to disorders of B-cell differentiation, proliferation, and survival,32–34 such as the NF-κB signaling pathways (PIM1, MYD88, TBL1XR1, PRDM1), epigenetic (KMT2D, KMT2C), cell cycle regulation (BTG2), MAPK/ERK signaling pathways (DUSP2), and NOTCH1.35–40 The NF-κB signaling pathway is a multi-component pathway regulating the expression of various genes, and the mutation of NF-κB signaling related genes was frequently observed in DLBCL.41,42 Interestingly, we found PIM1 was mutated in most PT-DLBCL patients (74%, Figure 3), while the mutation of NOTCH1 was mutually exclusive to PIM1 in our data, suggesting the variants of NOTCH1 may be another mechanism or a bypass event involved in PT-DLBCL.17,43 In addition, the MYD88 gene were mutated in nearly 50% of patients, whereas the concordant mutations of CD79B genes were not observed in our data, which may be caused by the lack of CD79B probe in exome sequencing panel.44 Notably, the BTG2 was frequently mutated in PT-DLBCL patients as well, which was consistent with previous research that missense mutations of BTG2 were relatively common in DLBCL,45,46 especially in patients with relapsed/refractory DLBCL (>20%).47,48

Clinical stage and IPI scores are classical indicators for prognosis estimation in DLBCL patients,49,50 and the gene mutations have been increasingly recognized to impact prognosis, eg, RAS‐MAPK pathway genes.51,52 We assessed the impacts of gene mutations in PT-DLBCL, as well as the pairwise mutation patterns. Notably, the mutations of BTG2 mutation was an independent prognostic factor in PT-DLBCL (poor OS prognosis, HR=2.23, 95% CI 1.15 to 4.33, P=0.007), while the mutation of TBL1XR1 and PRDM1 genes were risk factors as well (Table 3). Meanwhile, we found that the mutation of BTG2 could act as an important add-on for the prognosis of PT-DLBCL, even for patients with similar IPI score (Figure 7B). BTG2 as a transcription coregulator could meditate the activity of transcription factors,48 and serve as tumor suppressors in the MCD subtype of DLBCL by suppression of BLIMP1 (PRDM1).35,36 Furthermore, we investigated the expression profiles of BTG2 in TCGA DLBCL cohort, and found the expression of BTG2 gene was up-regulated in the comparison of tumor-vs-normal (Figure 5A), implying BTG2 was important molecule in tumor environment. Thus, the mutation of BTG2 may result in the dysfunction of tumor-repressor role of BTG2, and thereby contribute to malignant proliferation of B-cells in PT-DLBCL. In addition, we observed that 24% of patients were mutated in PRDM1 and were associated with shorter OS (P=0.040, Table 3). Besides, TBL1XR1 showed a significant co-mutation profile with BTG2 in our data (Figure 4), which may contribute to the tumor invasion in PT-DLBCL as well.39

Although this study had revealed the mutation profile of 360 lymphoma-related genes in 76 patients with PT-DLBCL, our data involved a limited number of different BTG2 mutations and only a small number of subjects with each specific mutation. More patients per mutation and a greater spectrum of BTG2 mutations are needed to be validate from multiple multi-center PT-DLBCL cohort and examined the clinical risks associated with mutations and immunohistochemistry for the BTG2. Thus, we had limited power to show the risk heterogeneity of patients with different mutation types of BTG2 and needed to be validated by multi-center data.

To the best of our knowledge, due to the rare subtype of DLBCL as PT-DLBCL, systematic analysis on the mutation map of gene was still lack in PT-DLBCL. Our research showed mutation profiles of lymphoma-related genes and their influence on prognosis of patients with PT-DLBCL mutations as well as several factors including age and IPI score etc. Our data demonstrated that the mutation of BTG2 could serve as a valuable prognosis indicator and candidate therapeutic target in PT-DLBCL.