Dissection of DLBCL microenvironment provides a gene expression-based predictor of survival applicable to formalin-fixed paraffin-embedded tissue.

Background: Gene expression profiling (GEP) studies recognized a prognostic role for tumor microenvironment (TME) in diffuse large B-cell lymphoma (DLBCL), but the routinely adoption of prognostic stromal signatures remains limited. Patients and methods: Here, we applied the computational method CIBERSORT to generate a 1028-gene matrix incorporating signatures of 17 immune and stromal cytotypes. Then, we carried out a deconvolution on publicly available GEP data of 482 untreated DLBCLs to reveal associations between clinical outcomes and proportions of putative tumor-infiltrating cell types. Forty-five genes related to peculiar prognostic cytotypes were selected and their expression digitally quantified by NanoString technology on a validation set of 175 formalin-fixed, paraffin-embedded DLBCLs from two randomized trials. Data from an unsupervised clustering analysis were used to build a model of clustering assignment, whose prognostic value was also assessed on an independent cohort of 40 cases. All tissue samples consisted of pretreatment biopsies of advanced-stage DLBCLs treated by comparable R-CHOP/R-CHOP-like regimens.
Results: In silico analysis demonstrated that higher proportion of myofibroblasts (MFs), dendritic cells, and CD4+ T cells correlated with better outcomes and the expression of genes in our panel is associated with a risk of overall and progression-free survival. In a multivariate Cox model, the microenvironment genes retained high prognostic performance independently of the cell-of-origin (COO), and integration of the two prognosticators (COO + TME) improved survival prediction in both validation set and independent cohort. Moreover, the major contribution of MF-related genes to the panel and Gene Set Enrichment Analysis suggested a strong influence of extracellular matrix determinants in DLBCL biology. Conclusions: Our study identified new prognostic categories of DLBCL, providing an easy-to-apply gene panel that powerfully predicts patients' survival. Moreover, owing to its relationship with specific stromal and immune components, the panel may acquire a predictive relevance in clinical trials exploring new drugs with known impact on TME.

Key Message

A 45-microenvironment-gene panel, derived from CIBERSORT deconvolution, effectively predicts survival of DLBCL patients using formalin-fixed, paraffin-embedded biopsy tissue samples. Its integration with cell-of-origin assignment improves prognostication.

Introduction

Reliable prognostic stratification at diagnosis and improvement of first-line/salvage strategies represent unmet clinical needs in diffuse large B-cell lymphoma (DLBCL). Gene expression profiling (GEP) studies suggested that patients’ outcome is entwined with a remarkable heterogeneity of both malignant clone and cellular/extracellular microenvironment, and unveiled prognostic cell-of-origin (COO) subtypes of DLBCL, namely GCB and ABC, driven by peculiar oncogenic pathways [1]. This prompted the adoption of immunohistochemical algorithms or NanoString-based assays, as GEP surrogates, for prognostic purpose [2].

Parallel studies emphasized the prognostic role of genes related to tumor microenvironment (TME) [3]. Lenz et al. demonstrated that the ‘Stromal-1’ and ‘Stromal-2’ signatures correlate with good and poor outcome in R-CHOP-treated patients, respectively, independently of COO [4]. However, these signatures failed in recognizing definite TME cytotypes with prognostic significance, and their practical use was limited by the lack of standardized assays applicable to formalin-fixed, paraffin-embedded (FFPE) samples. Despite several attempts [5-9], neither specific, nor reproducible TME biomarkers have been thus far identified, capable of predicting patients’ outcome.

Computational methods of GEP deconvolution allow high-sensitivity discrimination of cell subsets within complex tissues, as tumors [10]. Such approaches provide quantitative/functional information also on rare tumor-infiltrating elements, offering the unprecedented opportunity of reanalyzing available genomic data to explore correlation between TME cells and clinical outcomes.

Here, we applied the computational algorithm CIBERSORT to GEP datasets from pretreatment DLBCL biopsies to draw a map of immune and stromal components of TME. We identified microenvironmental prognostic genes and used the NanoString technology to develop a reproducible assay, which was validated in independent DLBCL cohorts.

Materials and methods

Patient cohorts

We selected two homogeneous cohorts of 175 and 40 patients with newly diagnosed, nodal DLBCL, not otherwise specified, from two multicenter trials (RHDS0305 [11] and DLCL04 [12]) and three monocenter ‘real-life’ selections from IRCCS - Istituto Tumori ‘Giovanni Paolo II’ of Bari (Italy), S. Orsola-Malpighi Policlinic, Bologna (Italy) and Hospital Clinic, Barcelona (Spain), respectively. Characteristics of the patients and selection criteria are shown in Figure 1A and detailed in supplementary methods, available at Annals of Oncology online.

Figure 1.

Characteristics of DLBCL patients in the study cohorts and overall study design. (A) Clinical data of patients form multicentric trials and ‘real-life’ cohort. (B) GEP data from 482 fresh-frozen DLBCL biopsies were analyzed by CIBERSORT to obtain a set of prognostic cytotype-related genes that were incorporated in a definitive 45-gene TME panel. The prognostic power of the panel was first assessed on 218 selected cases from the in silico cohort. Subsequent validation was carried out by NanoString technology on 175, FFPE samples from clinical trial and 40 ‘real-life’ patients. In addition, 79 randomly selected cases from the clinical-trial validation cohort were analyzed on a second NanoString Platform, as technical replicate. Based on the expression matrix from the 175 validation cases, a Random Forest classifier was built to assign each of the 40 ‘real life’ cases to a certain gene expression cluster and perform survival analysis. Thus, a composite model of survival prediction was developed by integrating the prognostic contribution of both TME and COO. FFPE, formalin-fixed paraffin embedded; CHT, chemotherapy; ASCT, autologous stem cell transplantation; COO, cell-of-origin; IPI, international prognostic index; TME, tumor microenvironment.

CIBERSORT analyses, NanoString-based gene quantification and building of a Random Forest-based model for survival prediction

A CIBERSORT-based deconvolution of GEP datasets (GSE10846 and GSE34171) from 482 DLBCLs was carried out using a 1028-gene signature matrix customized by normalizing referenced microarray data (Affymetrix) from 13 immune and 4 stromal cell types, according to CIBERSORT instructions (https://cibersort.stanford.edu/; supplementary Tables S1 and S2, available at Annals of Oncology online). DLBCL cases were stratified in ‘good’ and ‘poor’ outcome subgroups [according to overall survival (OS)] and only those cell types showing significantly different infiltration percentages between groups were analyzed. Then, the most expressed genes (log2 transformed >12) as well as those identified by a Random Forest analysis as the ‘top-30’ from each cytotype provided a selection of 150 genes (supplementary Figures S1 and S2, available at Annals of Oncology online). Among these, only genes showing a significantly different expression between ‘poor’ and ‘good’ subgroups were incorporated in a definitive prognostic panel of 45 TME-related genes (supplementary Figure S3 and Table S3, available at Annals of Oncology online). Further details are provided in supplementary methods, available at Annals of Oncology online. The prognostic performance of the panel was firstly assessed by an unsupervised clustering and survival analysis of 218 cases homogeneously selected (stages III–IV) from the GSE10846 dataset [4] (fresh/frozen biopsies). Gene expression levels were extracted from a normalized matrix and a long-rank test used to compare OS in the clustered samples.

Expression of TME genes and COO were assessed by the nCounter Analysis System (NanoString Technology) on the 175 FFPE cases enrolled in trials. Also, 79 of these cases were randomly selected as technical replicate using a second nCounter Analysis System located at a different Institute. Details are described in supplementary methods, available at Annals of Oncology online.

A Random Forest classifier was built on the expression of TME genes from the validation set and applied to 40 ‘real life’ cases to perform a TME clustering and survival analysis. Finally, a composite model of survival prediction was developed by integrating the prognostic contribution of both COO and TME data, and stratifying cases (both validation and ‘real life’ cohorts) into high, intermediate and low risk categories. The overall study design is outlined in Figure 1B.

Gene set enrichment analysis and in situ staining

Gene set enrichment analysis (GSEA) [13] was run on three independent GEP datasets (GSE10846 [4], GSE34171 [14], and GSE12195 [15]). Immunohistochemical and immunofluorescence studies were carried out on representative DLBCL cases (supplementary methods, available at Annals of Oncology online).

Statistical analysis and risk categories

Principal component analysis (PCA) was carried out to distinguish cell subsets in the customized matrix. Comparison between groups was carried out by independent t-test (two tails, unequal variance) and Mann–Whitney nonparametric test. PCA, Heatmaps, Kaplan–Meier estimator of survival and P values were produced using ‘R’ statistical software. Long-rank test was used to compare OS and progression-free survival (PFS) among patients in different groups. Analysis of parameters included distance calculations by the Euclidian methodology; Ward’s method was applied for cluster aggregation. Multivariate and univariate analyses were constructed through the Cox proportional hazards regression model.

Results

CIBERSORT-based enumeration of DLBCL-infiltrating cell fractions has prognostic significance

PCA confirmed that the customized matrix clearly distinguishes two main cellular groups including immune cells of both lymphoid and myeloid derivation, and cells of mesenchymal origin, respectively (supplementary Figure S1A, available at Annals of Oncology online). However, among the lymphoid component a variable degree of transcriptional overlap was expected between tumor and B-lineage cells, as previously described [16]. In fact, our analysis revealed that, among immune subsets, memory B cells (as expected in a B-cell malignancy), natural killer (NK) activated cells, dendritic cells (DCs), CD4+ T cells, and M2-type macrophages were the most represented cell fractions within DLBCL microenvironment, whereas myofibroblasts (MFs) prevailed among the mesenchymal cytotypes (supplementary Figure S1B, available at Annals of Oncology online). When correlated to clinical outcome, the cases with higher proportions of MFs, DCs, and CD4+ T cells showed significantly longer OS (P = 0.002, P = 0.02, and P = 0.004, respectively), whereas activated NK and plasma cells (PCs) were more represented in patients with poorer outcome (P = 0.001 for NK). Owing to the potential overlap between transcriptional programs of PCs and malignant B cells of DLBCL, PCs were not included in subsequent in silico analyses. Stratification according to median percentage of infiltration of MFs, DCs, and CD4+ T cells confirmed their in silico prognostic significance, since Kaplan–Meier curves indicated their direct correlation with a longer OS. On the contrary, activated NK cells correlated with a significantly shorter OS (supplementary Figure S1C, available at Annals of Oncology online).

Development of a 45-gene microenvironment panel

DLBCL cases were dichotomized according to the median infiltration percentage of cell types endowed with prognostic significance in silico. Subsequent analyses provided a selection of 30 MF-, 10 DC-, and 5 CD4+ T-cell-related transcripts, which was incorporated in a 45-gene microenvironment panel for the validation phase (supplementary Figures S2 and S3 and Tables S1–S3, available at Annals of Oncology online). Moreover, a GSEA showed a significant enrichment of genes biologically linked with cytotypes consistent with the high expression of genes in our MF-, DC- and CD4+ T-cell signatures (supplementary Table S4, available at Annals of Oncology online), demonstrating that they were representative of these prognostically relevant cell types.

The 45-gene microenvironment panel stratifies DLBCL patients on both fresh/frozen- and FFPE-derived RNA

Based on the expression of the 45 genes of the TME panel, an unsupervised clustering analysis efficiently stratified cases from the GSE10846 dataset in three distinct subgroups showing significantly different OS (P < 0.0001) (supplementary Figure S4, available at Annals of Oncology online). Using the NanoString technology on 175 FFPE biopsy samples of the validation cohort, the unsupervised hierarchical clustering analysis stratified cases into high, intermediate, and low gene expression clusters (Figure 2A) comprising patients with significantly different OS (Figure 2B) and PFS (Figure 2C). In particular, the cases with lower gene expression (cluster 3) showed significantly worse survival (median 5-year OS and PFS of 60.7% and 58.9%, respectively) than those at intermediate (cluster 2, 91.5% and 84.7%) and higher expression (cluster 1, 83.2% and 73.2%). When stratified according to the expression of specific cytotype-related genes, the unsupervised clustering generated subgroups with similar prognostic trend (supplementary Figure S5, available at Annals of Oncology online). In fact, Kaplan–Meier curves showed significantly different OS and PFS among expression clusters of MF-related genes (supplementary Figure S5A, available at Annals of Oncology online), whereas DC and CD4+ T-cell signatures identified subgroups with only significantly different OS (supplementary Figure S5B–C, available at Annals of Oncology online). Notably, a Cox multivariate model of OS indicated that the prognostic performance of the gene panel is independent of COO classification and international prognostic index (IPI) score (Figure 3A and supplementary Table S5, available at Annals of Oncology online). Finally, the analytical validity of these results was also confirmed by their inter-laboratory reproducibility supported by the analysis of 79 technical replicates using a NanoString platform placed in a different laboratory (data not shown).

Figure 2.

Prognostic subgroups of DLBCL based on microenvironment gene expression. (A) The heatmap depicts the unsupervised hierarchical clustering of 175 DLBCL cases (NanoString technology) and identifies three different clusters according to high (cluster 1), intermediate (cluster 2), and low expression (cluster 3) of all genes in the TME panel. The relative levels of transcripts are indicated according to the color scale. Each row group comprises genes associated to specific tumor-infiltrating cell populations and each column a biopsy sample. Kaplan–Meier curves of OS (B) and PFS (C) show that patients in clusters 1 and 2 have significantly longer OS and PFS than those in cluster 3.

Figure 3.

COO/TME combined model for survival prediction in DLBCL. (A) Forest plot of multivariable hazard ratios for OS. Multivariable analysis was adjusted for COO, clusters of microenvironment gene expression and IPI. (B) Kaplan–Meier curves showing OS (left panel) and PFS (right panel) of patients from clinical trial and (C) ‘real-life’ validation cohorts, according to the composite (COO/TME) survival risk model. Patients are assigned to one of the risk subgroups in dependence of their COO and TME categorization, as reported in the exemplificative color panel. High-risk category includes ABC patients belonging to cluster 3; intermediate-risk category comprises ABC patients assigned to cluster 1 or 2, and GCB or unclassified to cluster 3; and low-risk category contains GCB or unclassified cases belonging to cluster 1 or 2. N, number of patient; COO, cell-of-origin; IPI, international prognostic index.

Integration of TME- and COO-based data from 175 patients in the validation cohort (supplementary Figure S6, available at Annals of Oncology online) produced a survival model identifying three different risk categories. The high-risk category includes ABC cases belonging to cluster 3, showing the worst outcomes; the intermediate-risk category comprises ABC patients assigned to cluster 1 or 2, and GCB or unclassified to cluster 3, with longer OS and PFS; and the low-risk category containing GCB or unclassified cases belonging to cluster 1 or 2 (Figure 3B). The reliability of this model was also assessed on 40 ‘real-life’ cases assigned to a certain TME cluster by a Random-Forest-based model built on the validation set. Interestingly, the risk categories produced an efficient survival prediction also for real-life cases (Figure 3C), which expectedly showed a global worse outcome compared with clinical trial cases, probably since they lacked the same selection procedure.

Taken together, these findings indicate that the microenvironment gene panel retains an independent prognostic power that, in combination with COO, acquires a remarkable capability of survival prediction, identifying patients at extreme risk within either the GCB or the unclassified or the ABC subsets.

Potential microenvironment determinants of DLBCL prognostic categories

By applying GSEA, we found a significant enrichment of the ‘Stromal-1’ signature [4] in cases with higher MF infiltration [normalized enrichment score (NES) = 1.60, false discovery rate (FDR) = 0.06]. Such cases were also enriched in distinctive transcriptional programs of fibroblastic reticular cells (FRCs) and alpha-Smooth Muscle Actin (SMA)-expressing elements (NES = 1.88, FDR = 0.01 and NES = 2.5, FDR < 10−4, respectively) [17]. Consistently, the cases with higher alpha-SMA expression showed a significant enrichment of FRC genes (supplementary Figure S7A, available at Annals of Oncology online), and those with the highest MF infiltration were enriched in genes implicated in extracellular matrix (ECM) remodeling (supplementary Figure S7B, available at Annals of Oncology online).

We in situ analyzed the expression of proteins encoded by four fronting genes of the MF signature, namely Fibronectin, Collagen-I, Laminin, and Collagen-IV. Double-marker immunofluorescence on representative DLBCLs belonging to different clusters showed high variability in both distribution and intensity of ECM markers (supplementary Figure S8, available at Annals of Oncology online). In some cases, Fibronectin and Collagen-I formed a diffuse interstitial meshwork, while in others they mainly co-localized at perivascular foci. On the contrary, Laminin and Collagen-IV distribution varied according to the microvascular density. Also, the alpha-SMA staining revealed higher densities of MF infiltration in samples from cases with favorable outcome (supplementary Figure S8, available at Annals of Oncology online). Overall, the extreme variability of ECM protein expression by in situ immunostaining does not support its use as a reliable assay to surrogate the expression of the relative TME prognostic genes.

Discussion

In a quite dogmatic view of DLBCL composition, the paucity of TME suggests that malignant clones grow independently from microenvironment stimuli [18]. Only few GEP studies [4, 5] indicate that stromal signatures can predict survival, but the clinical application of such signatures is restrained by the lack of robust, reproducible and cost-effective biomarkers. In addition, the above-mentioned studies were mostly carried out on fresh or frozen samples, which are not available in most patients.

Herein, we applied a computational strategy to reanalyse GEP from whole-tissue biopsies and resolve in silico the cellular components of DLBCL. We found that a prevalence of MFs, DCs, and CD4 T cells and their gene signatures predict significantly longer survival after chemoimmunotherapy. From a clinical point of view, we capitalized GEP data to validate a new risk model that purely reflects the cellular composition of TME and categorize DLBCL patients in new COO-independent prognostic subgroups. Neither early-stage, nor high-grade ‘double hit’ patients have been included in our validation cohorts to avoid potential prognostic biases relying on low tumor burden or molecular aggressiveness of the disease. This emphasizes that current COO classification may fail in discriminating subsets of patients whose prognosis could be better predicted by a combined COO/microenvironment assessment. We validated a NanoString-based assay that provides simultaneous COO and TME information, and a Random-Forest model to facilitate the assignment of cases to a certain TME gene expression cluster in the clinical practice.

From a biological point of view, our evidence raises key questions on the nature and role of stromal elements in the pathophysiology of lymph nodes [19]. It is conceivable that cells recognized by CIBERSORT as MFs belong to a resident pool of FRCs, as immunologically specialized stromal elements [20]. Consistently, our signature includes genes encoding ECM molecules, as collagens (COL-1, -3, -4, -5, -6), proteoglycans and glycoproteins (Lumican, Byglican, Laminin, Fibronectin, SPARC), matrix metalloproteinases (MMPs), metallopeptidase inhibitors (TIMP1, TIMP3), adhesion molecules (ALCAM, AMICA1, Periostin), and regulators (TGFBI). GSEA indicated a significant enrichment of FRC-related genes involved in ECM organization in highly MF-infiltrated cases. The ‘Stromal-1’ signature [4] was also found enriched in these cases, suggesting that a specific pool of stromal cells, not clearly recognizable by GEP, may functionally influence the DLBCL biology. This is consistent with the attitude of stromal elements to react to inflammatory stimuli and remodel the microarchitecture of lymph nodes to orchestrate the immune response [17, 21]. Monocyte/macrophage precursors, B and T lymphocytes as well as other antigen-presenting cells (APCs), including DCs, are recruited within tumor milieu by gradients of chemokines, cytokines and growth factors [22] and traffic within the lymph node along collagen scaffolds mainly regulated by stromal cells [23, 24]. In turn, fibroblastic and immune cells participate to the local synthesis of soluble mediators active on APCs, effector cells and tumor cells themselves [25]. Some of the CD4+ T-cell- and DC-related genes in our panel, namely ALCAM, AMICA1, COL4A2, TNFSF13B, and SMAD1, encode activators of T lymphocytes and APCs via adhesion and paracrine mechanisms. Our findings suggest a functional relationship between MFs, ECM and immune system components in DLBCL, where the TME heterogeneity reflects unique paracrine and spatially organized interactions between tumor, mesenchymal and specialized subsets of immune cells [23]. This crosstalk, which may involve other lymphoid and myeloid cells as PCs, CD8+ T cells, NKs, and T-reg contributes to the immune control of the malignant clone and impacts the efficacy of chemo-immunotherapy [26, 27].

At a clinical level, the prognostic impact of CD4+ T cells and DCs in DLBCL has been previously reported [7], emphasizing the idea that patients with consensual activation of stromal and immune elements benefit from an empowered immune-surveillance in terms of survival. This study, however, remains uninformative about whether peculiar quantitative/functional assets of nontumor cells in lymph nodes precedes or follows tumor initiation and growth. Also, the role of some ECM proteins, such as SPARC, remains ambiguous since it may exert opposite influence on the tumor clone [28]. It controls the inflammatory milieu by remodeling ECM and its experimental deficiency promotes lymphomagenesis in vivo [29], suggesting again that peculiar stromal/immune programs may affect tumor behavior. Limitations of our study also include its retrospective nature and the small ‘real-life’ sample size. Further validation of the proposed model is required on large-scale, prospective cohorts of patients at different disease stages or selected according to the more recent molecular variants of DLBCL [30].

On the other hand, our findings indicate that the assessment of certain transcriptional programs, as those detected by our gene panel, may provide more reliable prognostic information than immunohistochemistry surrogates. We believe that computational approaches decoding ‘big’ genomic data, followed by proper validation, may capture functional processes linked to definite cytotypes within TME, and NanoString technology for gene quantification from FFPE tissues is particularly advantageous in term of feasibility, reproducibility and costs, and worthy of being routinely included in the diagnostic/prognostic workflow. Finally, although this work does not provide definite therapeutic targets related to TME, it prompts future investigations to assess whether our model can be predictive for patients treated by novel drugs with broad modulatory properties on TME.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.