Gene expression predicts overall survival in paraffin-embedded tissues of diffuse large B-cell lymphoma treated with R-CHOP

Diffuse large B cell lymphoma (DLBCL) is a heterogeneous disease with variable patient survival. It accounts for nearly 35% of all cases of lymphoma. Gene expression profiling (GEP) studies of DLBCL have been performed by different research groups and have identified largely nonoverlapping gene sets associated with patient survival. The Leukemia and Lymphoma Molecular Profiling Project reported on 17 genes that could be used to determine an outcome predictor score using a competitive microarray platform, that when divided into quartiles was able to predict patient overall survival.¹  A reanalysis of these data identified a “redox” signature score, which could also predict survival.²  Another group of investigators using a different microarray technology platform identified 13 different genes predictive of overall survival.³  A third group evaluated genes reported to be of prognostic interest in the literature to create a 6-gene model using quantitative reverse-transcribed polymerase chain reaction (RT-PCR).⁴  The genes identified by these different researchers represent largely nonoverlapping gene sets with the exception of 2 genes, BCL6 and fibronectin, as having prognostic significance. No actual side-by-side comparisons have yet been published. These conflicting data in the literature make it difficult to determine which genes are the most prognostically important to evaluate in new studies. Furthermore, progress is often limited by the numbers of cases for which frozen materials are available.

To overcome the limitation of using snap-frozen tissues, we used a multiplexed quantitative nuclease protection assay, the ArrayPlate, qNPA, useful for measuring mRNA levels in fixed paraffin-embedded samples. The assay was customized to measure all of the genes of interest from 4 previously described gene expression papers of DLBCL.⁵  This assay's performance demonstrated excellent reproducibility, applicability to archived paraffin blocks, and quantitative results that correlated well with GEP.

Most patients are now treated with monoclonal antibody therapy, most commonly rituximab, combined with chemotherapy, raising the question of whether the prognostic genes identified in the setting of chemotherapy treatment alone retain prognostic significance. The results of several randomized trials and a population-based registry experience have clearly indicated that rituximab plus cyclophosphamide, hydroxydaunorubicin, oncovorin (vincristine) and prednisione (R-CHOP) is the new treatment standard for DLBCL.^6-8  There is some evidence to indicate that, indeed, the importance of some factors may be affected with new treatment, in particular BCL6 and BCL2. In one study, the authors found a reduction in treatment failure and death with the addition of rituximab to CHOP only for the BCL6⁻ cases and that addition of rituximab did not benefit BCL6⁺ cases.⁹  Another group demonstrated that the addition of rituximab to chemotherapy for DLBCL patients overcame the negative prognostic value of the BCL2 protein.¹⁰  Thus, the question of how combined immunochemotherapy modifies the prognostic ranking of specific genes measured by expression levels remains pertinent.

In this paper, we demonstrate the robust use of the qNPA assay on formalin-fixed, paraffin-embedded tissue (FFPET) blocks, that the results can be related to patient overall survival, that prognostic genes previously identified remain relevant in the R-CHOP era, and that loss of immunosurveillance and high proliferation together identify patients with the worst outcome.

Three- to 5-micron unstained cuts from FFPET blocks were used from 93 cases of DLBCLs treated primarily with CHOP or similar CHOP-like chemotherapy and 116 cases treated with R-CHOP. Cases of transformed lymphomas were excluded. Frozen blocks from the CHOP-alone cases had been analyzed as part of a prior publication.¹  As previously reported, these cases had undergone consensus review by a panel of expert hematopathologists and confirmed as DLBCL. The R-CHOP cases were taken from the current case files at the University of Arizona, British Columbia Cancer Agency, and Oregon Health Sciences Center. Of these 116 R-CHOP cases, frozen blocks from 32 were also used in another study and had undergone review by an expert panel.¹¹  All tissues used for this retrospective study came from pretreatment diagnostic biopsies using excess diagnostic tissue under institutional review board–approved protocols. This project has been reviewed and approved by the University of Arizona Human Subjects Institutional Review Board, in accordance with the Declaration of Helsinki.

The performance of the ArrayPlate assay customized for use in DLBCLs has been described previously by our group.⁵  Additional descriptive details of the assay, schematic figure, and probe sequences are given in Table S1, Figure S1, and Document S1 (available on the Blood website; see the Supplemental Materials link at the top of the online article). The complete previous publication is available at the online archives for Laboratory Investigation.⁵ 

We used the key genes¹²  identified as prognostically important in 4 previous papers of DLBCLs, which accounted for 36 genes of interest,^1-4  as listed in Table 1 in the order in which they were listed in the original references. The housekeeping gene, TATA Box Binding Protein, was chosen for normalizing the data based on its stable expression at moderate levels in 12 lymphoma cell lines and 80 B- and T-cell lymphoma samples compared with 11 other “housekeeping” genes, using quantitative RT-PCR as well as our own previous experience with this gene in the ArrayPlate assay, which showed it to be moderately expressed with minimal variability in all samples tested to date.^5,13 

Statistical analyses of the association of gene expression, as measured by Array Plate qNPA technology, with survival were performed on the 116 cases treated with R-CHOP and the 93 cases treated with CHOP or CHOP-like regimens alone. The logarithms of gene-expression values were standardized to have standard deviation equal to 1.

Hazard ratios, 95% confidence intervals, and P values for the univariate associations between standardized log gene expression levels and patient overall survival were obtained using Cox proportional hazard regression.¹⁴  To account for the relatively large number of test statistics, the overall statistical significance of the set of hypothesis tests against a global null hypothesis of no association was calculated, by permutation resampling, based on the “tail strength” (TS) statistic.¹⁵  A test of the overall statistical interaction of gene expression and treatment type was considered in a similar fashion.

In an exploratory analysis of parsimonious multivariate models, we used best subset selection, which determines the “best” model based on the global score χ² statistic.¹⁶  Candidate genes used in the model building process were those achieving nominal P values less than .05. The top 3 models for each of 1, 2, 3, and 4 variable models were derived. For presentation purposes, for each factor included, the overall best identified model, patients were categorized by high versus low gene expression (above or below the median value). Patients were then grouped according to the number of adverse risk factors, and survival was examined.

Finally, we assessed whether the AP gene-risk model retained significance of the biologic aspects of the malignant cells after adjusting for the clinically based IPI index.¹⁷ 

Separately, cutoff point analysis was performed on the factors identified in multivariate modeling to optimize identification of expression levels of highest risk. Permutation resampling was used to adjust significance levels of the proportional hazards score tests among all evaluated cutoff points.¹⁸  In addition, to control statistical variability of the cutoff point analysis, a minimum possible group size of 10% of total patients was set for the analysis.

Of 209 cases attempted, there was only one that did not result in adequately detectable signal. In situ hybridization using a polyDT probe (Ventana Medical Systems, Tucson, AZ) demonstrated that the mRNA was degraded (data not shown). This failure was therefore attributed to sample inadequacy rather than a technical failure of the assay itself. TATA-binding protein was moderately and consistently expressed in all samples, as in our previous work, and again used as the control gene for normalization of the data.⁵ 

Initial evaluation of HTG results included univariate analysis of individual gene levels with respect to patient survival in both treatment groups, using the logarithm of the gene expression measurements. To further explore potentially important genes, we calculated hazard ratios of death and assessed whether they trended in the direction predicted by the previously reported literature. For each of the treatment groups, we assessed the overall significance of the panel of genes using the tail strength statistic and permutation resampling. Any gene that was significantly associated (P < .05) with overall survival in univariate modeling was assessed for potential inclusion in a multivariate risk model using Cox regression analysis.¹⁴  To determine the best model, we used best subset selection, which determines the “best” model based on the global score χ² statistic, and adjusted this model for clinical IPI score.¹⁶  Lastly, a variable cutoff point analysis was performed on the 2 key genes to determine whether there were more relevant cutoff points, rather than the preselected 50th percentile, which might have biologic implications. Permutation sampling was used to adjust for multiple comparisons in the cutoff point optimization.

For chemotherapy-alone (mainly CHOP) treated cases, gene expression levels were significantly correlated with overall survival at P less than .05 for 15 of 36 prognostic genes, including the major histocompatibility complex (MHC) class II genes HLA-DRA and HLA-DPA1, germinal center–associated genes BCL6, GCET1 (SERPINA9), stromal associated genes (ACTN1, COL3A1, CTGF, FN1), proliferation genes MYC, CCND2, PRKCB1, as well as PDCD4, TLE, B4GALT1, and BCL-2. These genes represented all 4 prognostic signatures from Rosenwald et al,¹  4 of 13 genes reported by Shipp et al,³  and 3 of 6 genes from Lossos et al.⁴  An additional gene, CCL3, was borderline significant at a P value of .062.

For the R-CHOP–treated patients, 11 of the 36 genes analyzed were significantly associated with survival at the P less than .05 cutoff level. These genes were GCET1 (SERPINA9), HLA-DQA1, HLA-DRB, ACTN1, COL3A1, PLAU, MYC, BCL6, LMO2, PDCD4, and SOD2. An additional gene, FN1, was marginally significant at a P value of .078. Results of univariate analyses compared for the 2 treatment eras are shown side by side in Table 2. To emphasize the genes with recurrent significance, the P values at .05 or less are footnoted. Average 2-year overall survival for each gene cut at above and below the median expression level are also summarized in Table 2. We note that survival rates at 2 years were chosen as simple descriptive summary statistics. Similar results were seen with 3- and 4-year rates, although estimates were more unstable because of more censored cases. The P values presented in the tables are based on Cox score tests using the continuous logarithm of gene expression and therefore do not depend on the choice of summary survival estimates presented.

Comparative overall survival curves in the different treatment eras for HLA-DRB (an MHC class II gene), BCL6, and MYC are demonstrated in Figure 1. These examples demonstrate the ability of the ArrayPlate assay to generate meaningful quantitative data that can be related to patient outcome. The results also demonstrate that, for these well-known prognostic genes, there is continued evidence of prognostic relevance in R-CHOP–treated patients.

For most genes in both treatment groups, the estimated hazard ratios of death trended in the direction predicted by the original studies (Table 3). Hazard ratios correspond to a change in 1 standard deviation in log expression levels and a hazard ratio more than 1 indicates an association between high expression (above the median) with poorer outcome, whereas hazard ratios less than 1 indicate an association between high expression with better outcome. Therefore, an estimated hazard ratio that is very small in magnitude (eg, close to 0) corresponds to a gene with strong association between higher expression and longer survival.

To address the testing of the multiple genes in the panel, an overall test of the 36 P values was performed using the TS statistic and permutation resampling. There is evidence of association between the overall 36-gene panel and outcome in both CHOP-treated patients (TS, P = .007) and R-CHOP patients (TS, P = .013).¹⁵  We also considered an overall test of differences by treatment group in the association between each gene and survival (statistical interaction). Although power for interaction testing is limited, there was no evidence of a differential effect of the overall 36-gene expression panel between the 2 treatment types (TS, P = .250).

As an overall assessment of important prognostic features, the IPI distribution was assessed among patients in the 2 treatment types. In the CHOP-alone patients, 41% had IPI of 0 to 1, 48% had IPI of 2 to 3, and 11% had IPI of 4 to 5. In the CHOP-R-treated patients, 40% had IPI of 0 to 1, 56% had IPI of 2 to 3, and 4% had IPI of 4 to 5. There was no evidence of a difference between the 2 treatment groups (P = .18). Table 4 details the distribution of the individual factors of the IPI score between patients in the 2 treatment eras.

As an exploratory analysis of multivariable prognostic models, best 1, 2, 3, and 4 variable models, as determined by best subsets analysis, were calculated (data not shown). The best 2-variable model was the combination of MYC and HLA-DRB, with a model χ² of 16.6. However, we note that other 2-variable models, including MYC with HLA-DQA1 or PLAU, had modestly smaller model χ² statistics. Given the relatively small number of events in this study, conclusive statements about the overall best model are not possible. There was no evidence that 3 variable models yielded any statistical improvement in model fit.

Patients were defined as having high or low levels of MYC and HLA-DRB. Twenty-eight patients (24%) had both adverse gene levels. These patients had much worse survival than patients with 0 or 1 adverse gene level (2-year overall survival, 38% vs 87%), as shown in Figure 2A. Differences are presented for both the high and low IPI subgroups (Figures 2B,C). The survival disadvantage for patients with both adverse gene levels appears particularly pronounced in patients with high IPI (2-year estimate, 14% vs 68%), although there was no evidence of an interaction between number of adverse gene levels and IPI group (P = .88).

Both CHOP and R-CHOP data were combined to further explore the nature of the association of expression of these 2 genes with survival using cutoff point analysis. For HLA-DRB, the highest χ² value indicating the most significant cutoff point was at the 20th percentile of gene expression (P = .01 based on permutation resampling to account for the multiple testing). For MYC, the most significant cutoff point (P = .01) was at the 80th percentile of expression (Figure 3). Given the adaptive nature of cutoff point selection, any multivariate model based on cut-off point levels identified in this analysis would be best validated independently. We emphasize that, although the 80th percentile was the optimal cutoff point for MYC (corresponding to a χ² value of > 15 and a nominal P < .001), there were a wide range of cutoff point values that were also nominally significant (P < .025). This indicates that other cutoff points may lead to interesting prognostic models.

GEP studies in cancer have relied on snap-frozen tissues, which are often not obtained at diagnosis and are becoming decreasingly available in the current era of small biopsies. To overcome this obstacle, in this paper we describe the use of a quantitative S1 nuclease protection assay to measure gene expression levels in FFPET blocks of DLBCLs, related to results to patient outcome in 2 different treatment eras, and used multivariate and variable cutoff point analysis to identify molecules related to antigen presentation and cell cycle/proliferation as key indicators of patient outcome.

The ArrayPlate qNPA multiplexed mRNA assay used for this project was originally designed to support genomics-driven drug discovery efforts.¹⁹  The standardized protocol is automated and optimized for processing large numbers of samples. The S1 nuclease digestion that the qNPA incorporates is a well-established approach that predates more widely used PCR-based procedures for quantifying gene expression.²⁰  For a 12-gene signature predictive of phospholipidosis in human liver cells to examine treatments with 80 different compounds, mRNA scores correlated significantly (R² = 0.95) between the ArrayPlate and real-time PCR assays.²¹  Unlike RT-PCR, however, the ArrayPlate qNPA did not require RNA isolation, reverse transcription, or amplification. Instead, target mRNA-specific 50-mer qNPA probes were retained in a stoichiometric fashion. The direct probe-to-target hybridization on which the assay relies was critical for the successful GEP of FFPET blocks of material.⁵  We demonstrated that the assay measured mRNA, whether it was released from FFPET sections or cross-linked in situ. The use of relatively short target 50-mers in the sample meant that robust assay results were obtained even for archived specimens nearly 20 years old in which RNA degradation was presumably substantial, resulting in an extremely low failure rate of only 1 in 209 cases (0.5%). With most comparable techniques, failure rates of all attempted cases are seldom reported and anecdotally high.

In our series of CHOP-treated cases, qNPA results for many genes were significantly associated with overall survival when using a standard cutoff point at the median level of expression. These genes represented key expression signatures previously identified as prognostically important in large discovery-oriented GEP experiments. Because these cases had been a part of the original case series by Rosenwald et al, the results showing significant correlations in 10 of the 17 outcome predictor score genes from that paper served to demonstrate the ability of the qNPA to generate data from FFPET blocks that can be correlated to patient survival.¹  We also concluded that many of the genes identified by the different research groups each do represent important aspects of DLBCL biology, which deserve further exploration. We wish to emphasize that the biologic importance of these findings lies in the functional groups of genes rather than specific genes themselves because, in any one relatively small dataset, one or the other of a set of genes may reach statistical significance.

In the R-CHOP–treated patients, again many genes were significantly associated with overall survival using a cutoff at the median level of expression. For the genes that did not reach statistical significance, the trends in the hazard ratios largely agreed with those from the original publications indicating that many of the “nonsignificant” genes may become significant if analyzed in a larger dataset. Nonetheless, our findings indicate that representatives from each of the previous functional groups of genes remain relevant in the R-CHOP era and that the key biologic aspects of DLBCLs are not selectively affected by addition of monoclonal antibody to CHOP chemotherapy. This interpretation is supported by the global test of interaction between the 2 treatment groups, which did not show any significant difference in gene levels as related to outcome overall.

In agreement with our findings of the continued importance of expression levels of germinal center associated genes (BCL6, GCET1 in the CHOP era; and BCL6, GCET1, and LM02 in the R-CHOP era), a recent paper demonstrated that LM02 immunohistochemical staining was a significant predictor of survival in R-CHOP–treated patients similar to chemotherapy alone.²²  Recent reports using quantitative RT-PCR and GEP have described similar findings of the continued importance of most genes.^11,23  A possible exception is the BCL2 gene, which was significant in the CHOP era but not in our series of R-CHOP–treated cases (although the hazard ratio trended in the correct direction with higher BCL2 levels representing a higher risk). Whether this represents the real impact of adding rituximab to the chemotherapy regimen as reported by Mounier et al using immunohistochemistry needs to be clarified in a larger cohort of patients.¹⁰  Small changes in the significance of genes in CHOP- versus R-CHOP–treated patients would be difficult to find without large numbers of cases because most patients are simultaneously treated with rituximab and CHOP, which may obfuscate the survival changes from rituximab. Application of the variable cutoff point analysis (applied here to HLA-DRB and MYC only) to all of the genes in this study may be an interesting exercise to determine whether the use of different cutoff points allows the impact of gene expression to be seen more readily. However, we are not attempting to call out specific genes but rather to emphasize that the same groups of genes remain relevant in the current treatment era.

We and others have previously demonstrated that, with decreasing MHC class II expression, patient outcome worsens with those in the lowest 10% having an extremely poor outcome.^24-26  In this paper, we again demonstrated the significance of the MHC class II antigen gene expression levels (when cut at the median in our screening analysis). In addition, we show the strong significance of very low levels of MHC class II using the variable cutoff point analysis, which identified the most significant impact on survival at the 20% expression level consistent with our previous work. Although the MHC class II signature effect demonstrates similar trends with CHOP and R-CHOP, there is nonetheless a diminished MHC class II effect in the R-CHOP group, which could reflect some partial reversal of the negative prognostic value of MHC class II. A large study may be needed to further investigate this issue. The importance of lost MHC class II appears to be related to a loss of tumor immunosurveillance, as evidenced by decreased numbers of tumor-infiltrating T cells as reported by us and others.^24,27-29  Continued significance of these genes in the R-CHOP era is not surprising given the strong influence of tumor immunosurveillance on lymphomagenesis and outcome. This class of molecules therefore constitutes a therapeutic target in a subset of DLBCL.

Previous studies have identified high levels of MYC as being prognostically important in DLBCL in the CHOP era. It is therefore of interest that we were able to confirm the importance of MYC in the R-CHOP era and show, by using a quantitative method and variable cutoff point analysis, that the highest gene expression levels were most significant, and to narrow the collection of genes down to a useful signature of MYC and one other gene that identified a subgroup with only a 14% overall 2-year survival. This subgroup actually had approximately 40% of the observed deaths occurring within the first 2 years. MYC has a complicated story in aggressive B-cell lymphomas. In studies of DLBCLs, the presence of dual translocation of both BCL2, such as in t(14;18) as well as MYC as in t(8;14), have been identified as having poor outcome.³⁰  It is possible that dual translocations, alternative MYC translocations, amplifications, or mutations could all result in increased gene expression levels in these cases. Some of these cases could be DLBCLs that harbor a MYC translocation as a primary pathogenetic event, for which the impact on outcome is unknown. This is an excellent topic for future study.

The best selected 2-variable model predicting overall survival in DLBCL was the combination of MYC and HLA-DRB. It is striking that MYC-based proliferation and MHC antigen expression are important in R-CHOP–treated DLBCLs and are also the defining attributes of Burkitt's lymphoma suggesting that lack of immunosurveillance (through loss of either major histocompatibility class I or class II expression) and high proliferation are the key features defining the most aggressive B-cell lymphomas.^1,31,32  Down-regulation of immune response has been identified as a key feature of Burkitt lymphoma and poor-outcome DLBCL previously.^{24,28,29,31,33-36}  Highly proliferative DLBCL has also been identified as a key feature of poor-outcome DLBCL.^1,4,37  This tentative 2-gene model therefore is biologically plausible based on previous literature. Although these genes were identified using the median cutoff, the variable cutoff point analysis went on to demonstrate that those patients in either the lowest 20% of MHC class II expression or highest 20% of MYC expression have the poorest outcome. In addition, we note that the use of both high-risk IPI and the 2 genes identifies a subpopulation of 14 patients in which only 2 survive and 12 die during the first 2 years, yielding a 2-year overall survival rate of 14%. The remaining patients, with neither of these factors, appear to represent the curable population of DLBCL, whereas those with one or more of these factors should be the focus of future therapeutic development.

More broadly, the ability to use unstained cuts from fixed paraffin-embedded tissue blocks for multigene quantification has the potential to change the practice of pathology. In particular, the ability to interrogate routinely handled tissues, using minimal amounts of material, without RNA extraction, allows application of newly discovered gene expression biomarkers to patients uniformly treated on clinical trials (because fixed tissue is usually the only type of tissue easily available in such a context). Thus, this portable technology solves the major practical matter of obtaining and shipping snap-frozen tissue biopsies for analysis, allowing for the necessary validation studies required for biomarker application to clinical practice.

The authors thank the entire Leukemia and Lymphoma Molecular Profiling Project research group for producing the large well-annotated datasets that provided the data on the original group of CHOP-treated patients.

R.D.G. and J.M.C. are supported by the NCIC Terry Fox Program Project (grant 016003). L.R.'s laboratory is supported by the American Cancer Society (grant RSG0605501L1B).

Contribution: L.M.R. designed research, contributed cases, analyzed data, and wrote the paper; M.L.L. designed research, analyzed data, and wrote the paper; J.M.U. analyzed data and wrote the paper; T.P.M. and D.O.P. designed research and edited the paper; T.M.G. designed research and wrote parts of the paper; C.M.S. performed the research and analyzed data; R.R.M. contributed new analytical tools and edited the paper; B.S. supported new analytical tool development; R.M.B., E.C., A.R., L.H.S., J.M.C., and N.J. submitted case materials and reviewed the paper; and R.D.G. submitted cases, designed the research, and reviewed the paper.

Conflict-of-interest disclosure: R.R.M., C.M.S., and B.S. are employees of High Throughput Genomics. R.R.M. and B.S. hold stock and stock options in High Throughput Genomics. C.M.S. holds stock options in High Throughput Genomics. The remaining authors declare no competing financial interests.

Correspondence: Lisa M. Rimsza, Department of Pathology, University of Arizona, 1501 N. Campbell Avenue, PO Box 245043, Tucson, AZ 85724-5043; e-mail: lrimsza@email.arizona.edu.