A molecular risk score based on 4 functional pathways for advanced classical Hodgkin lymphoma

Classical Hodgkin lymphoma (cHL) is assumed to be a curable tumor, but an important fraction of patients with advanced disease do not respond favorably to the current standard chemotherapy regimens whose base is adriamycin. The most widely used and reproducible prognostic score is the product of clinical and analytical parameters integrated in the International Prognostic Score (IPS), but it still fails to identify accurately, at the moment of diagnosis, a significant fraction of patients with very poor prognosis.^1-3  Thus, the identification of biomarkers that, at diagnosis, may be consistently associated with nonresponse is essential for the recognition of patients at high risk of treatment failure to establish a more rational risk-adapted treatment strategy.

cHL represents a distinctive model of histologic complexity, with a minor population of the neoplastic Hodgkin and Reed-Sternberg (HRS) cells diluted in a reactive inflammatory background composed of nonneoplastic B and T cells, macrophages, eosinophils, neutrophils, and plasma cells. The complex relationship between the HRS cells and their microenvironment is only partially understood; however, important, if fragmentary, advances in our understanding are steadily being made.⁴  The clinical outcome of cHL has been found to be related to the expression of multiple biologic markers alone^5-8  or in combination,⁹  expressed either by the tumor HRS cells, macrophages, regulatory T cells, or other nonneoplastic cell subpopulations.^10-13 

Some of these previous analyses rely on array-based gene expression analyses, which use frozen tissue and in most cases can only provide retrospective information. Investigators of other studies have used immunohistochemical staining, with some inherent limitations to the reproducibility of the data thus generated. It is now feasible to apply multigenic predictive molecular tests in advanced cHL patients in a routine setting by the use of a quantitative reverse-transcription polymerase chain reaction (RT-PCR) assay as we here describe, incorporating a selected number of genes that capture information from tumor and microenvironment cell components, designed for application to routine formalin-fixed paraffin-embedded (FFPE) samples and that can be used at the moment of the initial diagnosis.

Previous studies allowed us to identify a group of genes whose expression was associated with the response of patients with advanced cHL to standard first-line treatment. Thus, the selection of genes to be analyzed was determined by results previously obtained in 2 independent series of 29 and 52 advanced cHL patients.^10,11 

The patients included in this study fulfilled stringent, previously described criteria (ie, age older than 16 years; advanced cHL; Ann Arbor stage IV, III, or IIB with bulky masses¹ ; proven HIV-negative status) who have been treated with a first-line standard chemotherapy regimen that included adriamycin—ABVD (adriamycin, bleomycin, vinblastine, and dacarbazine) or ABVD variants—and for whom information was available about the achievement of complete remission (CR) and a follow-up of at least 12 months thereafter, which is a well-known and accepted surrogate indication of the course of the disease. With respect to the latter, patients were considered to have had a favorable course if they had achieved CR and maintained it for at least 12 months or with an unfavorable course if they had either not achieved CR or if they had once had it but had relapsed during the following 12 months. All tissue samples consisted of representative pretreatment lymph node biopsies collected after revision and approved by the institutional review board of the participating institutions of the Spanish Hodgkin Lymphoma Study Group. The study initially included 282 FFPE patients, who were randomly split and assigned to the training (194 cases) or validation sets (88 cases) on the basis of the minimum estimated sample size to derive the final model (Table 1).

Additional exclusion criteria were insufficient RNA quality (purity ratio A260:A230 < 1.7) or a weak RT-PCR signal (average cycle threshold > 35) for the reference genes or in more than 10 genes of the assay. As a result, 20 patients were excluded, and the remaining 262 patients (183 in the training group and 79 in the validation group) meeting these criteria were included in the statistical analysis (Table 1).

The genes included in the assay initially were selected from 2 preliminary expression-profiling studies^10,11  that rendered a list of genes expressed by HRS and microenvironment cell subpopulations identified in unfavorable cHL patients. Selected genes were primarily chosen on the basis of their prognostic ability and capacity to represent biologic functions identified as relevant in cHL pathogenesis.¹¹  In addition, the strength and consistency of primer and probe performance also were taken into account.¹¹  The initial selection consisted of 30 genes, including genes expressed by the neoplastic cells involved in the cell cycle (G₂/M), apoptosis, histones, chaperones, drug metabolism, and mitogen-activated protein kinase signatures, and from microenvironment genes expressed by different cellular or functional populations of T cells, monocytes, macrophages, and dendritic cells. Details of the RT-PCR assays are available in supplemental Table 1 (available on the Blood Web site; see the Supplemental Materials link at the top of the online article).

Gene expression was analyzed by the use of a customized TaqMan low-density array platform (Micro Fluidic Cards; Applied Biosystems) on FFPE as previously described.^11,14  A preamplification step (PreAmp; Applied Biosystems) was used to improve the sensitivity of our assay for low-abundance target genes available from FFPE samples.^15-17  Reactions were performed by use of the ABI PRISM 7900HT Sequence Detection system (Applied Biosystems), and we measured the expression of each gene in triplicate and then normalized it with a set of 2 reference genes (HMBS and GUSB) whose uniform expression in cHL tumor samples was tested in previous studies.¹¹  Missing values were imputed using the K-nearest neighbor algorithm.¹⁸ 

Differences in the distributions of standard clinical parameters (age, sex, stage, IPS, the individual variables contained in IPS, and outcome) in the estimation and validation datasets were tested by the Pearson χ² test (Table 1).

The first end point of this study was the response to standard first-line treatment considering favorable response (F) and unfavorable response (U), as mentioned previously. Data from second-line and salvage therapies and/or bone-marrow transplantation were not considered.

The selection of the best predictive genes and the logistic regression model was on the basis only of the data from the training group of 183 patients, without any previous survival analysis that used information from the validation group. Univariate regression analysis was performed with treatment response (F vs U) as the dependent variable to identify genes significantly associated (P < .05) with outcome. In addition, final gene selection analysis was performed by cross-validation with the use of 3 prediction algorithms (http://tnasas.bioinfo.cnio.es/): diagonal linear discriminant analysis,¹⁹  support vector machines,²⁰  and K-nearest neighbor.²¹  Cross validation was used to test the classification ability of the initial set of significant genes to choose the strongest predictor genes, which were classified into functional groups on the basis of their known biologic relationship and their coregulated expression as estimated by the Pearson correlation coefficient. Individual genes from each functional group were weighted by the use of linear discriminant analysis.²²  Finally, these functional gene clusters associated with cHL outcome were analyzed in a multivariate logistic regression model with response to therapy (F vs U) as a dependent variable. In this way, an algorithm was derived that combines these measurements into a quantitative “molecular risk score” (MRS), which can be used as a continuous variable to estimate the probability of treatment response. The MRS cut-off points were prespecified by the use of an area under the receiver operating characteristic curve (ROC) analysis to define different risk groups. (See the supplemental Appendix for details of the statistical analysis and methods.) Finally, performance of the logistic regression model was tested in the validation group of patients (n = 79).

For graphical representation, survival analyses were performed with the Kaplan-Meier method and long-rank test separately in the training and validation series and in the entire series. Because the primary objective of the study was to identify patients at high risk of treatment failure, we used failure-free survival (FFS) as the fundamental end point for survival analysis. FFS was defined as the time interval between treatment initiation and treatment failure or last follow-up. Failure was defined as either the failure to achieve CR or the occurrence of progressive disease, irrespective of whether there had been an initial CR. Overall survival (OS), an end point whose significance is imperfect because it is conditioned by the effect of subsequent eventual treatments and complication of treatment, was included as a secondary end point in the survival analyses, defined as the time interval between diagnosis and death caused by the lymphoma.

Finally, in the whole series, a multivariate Cox proportional hazards model, including the data at diagnosis, the IPS stratified as previously defined (0-2 vs ≥ 3),¹  and its 7 individual variables (hemoglobin < 1.5 g/dL; albumin < 4 g/dL; leukocytosis ≥ 15 000/mm³; lymphopenia < 600/mm³; age ≥ 45 years; male sex; stage IV), was applied to test the independence of the MRS, including the remaining significant variables in a final integrative model.

All statistical analyses were 2-sided; values of P less than .05 were considered to be significant. These were performed with SPSS 15.0 (SPSS Inc). Survival curves were assessed by the Kaplan-Meier method, and risk groups were compared by the log-rank test. Plots were generated with the use of GraphPad Prism Version 5 (GraphPad Software, Inc).

In training series, univariate regression analysis of the expression data for the 30 initially selected genes revealed 20 genes to significantly predict failure to first-line treatment (supplemental Table 3). When cross validation was applied, the genes most frequently found in prognostic models consisted of a panel of 11 genes that were included in the final model: BCL2, BCL2L1, CASP3, HMMR, CENPF, CCNA2, CCNE2, CDC2, LYZ, STAT1, and IRF4.

To derive the model, we took a 2-step approach, first combining individual gene-expression patterns into precise functional pathways and then subsequently correlating these functional groups with the clinical outcome by using multivariate logistic regression. Final selected genes were weighted by the use of linear discriminant analysis and clustered into their corresponding functional pathways defined as macrophage activation (LYZ, STAT1), cell cycle (HMMR, CENPF, CCNA2, CCNE2, CDC2), and apoptosis (BCL2, BCL2L1, CASP3; Figure 1A). The Pearson correlation coefficient was significant for the genes included in each of the signatures (P < .001) apart from IRF4, which was included as an independent predictive gene because there were neither distinct functional relationships nor statistically significant correlations with other genes or pathways (supplemental Table 4). These functional groups captured information about the tumoral HRS cells and their nontumoral microenvironment, in agreement with previous studies.¹⁰  The multivariate logistic regression analysis integrating these pathways showed that cell cycle and apoptosis terms were associated with an unfavorable outcome of patients, whereas macrophage activation and IRF4 signatures had protective effects.

Thus, the optimized final model was determined on the basis of the relative contributions of each of the 4 functional terms, as described in the following equation: constant (−0.913) + (0.401 × apoptosis) + (0.284 × cell cycle) + (−0.301 × macrophage activation) + (−0.143 × IRF4). The continuous probability function generated by the logistic regression was defined as the MRS to treatment failure and ranged from 0.06 to 0.813 (Figure 1B). ROC analysis was used to define a threshold for stratifying patients, and the largest area under the curve was obtained by the use of 0.3 as the threshold, thus dividing the series into high-risk (MRS ≥ 0.3) and low-risk (MRS < 0.3) cases (Figure 1B-C).

The MRS differed significantly between the various outcome groups (supplemental Figure 2) and predicted treatment response with an accuracy of 68.9% in the estimation dataset and 70% in the validation dataset. Because treatment failure was the main end point used to derive the logistic regression function, FFS as previously defined, was used for validation of the model and Kaplan-Meier analysis of survival. Predicted probabilities also identified 2 risk groups associated with FFS in both the estimation and validation sets (Figure 1C; supplemental Figure 3). FFS probabilities at 5 years were 74.0% versus 44.1% (P < .001) in the training set and 67.5% versus 45.0% (P = .022) in the validation set.

In addition, the analyses of OS showed different risk groups identified by the MRS in the estimation group of patients. The differences were not significant in the validation series, probably because of the limited number of events (supplemental Figure 3).

There were no significant statistical differences in the MRS distributions of the 2 histologic subtypes, nodular sclerosis and mixed-cellularity cHL (P = .186, t test), and the differences in survival between risk groups identified by the MRS remained significant in both the estimation and validation sets stratified by histologic subtype (supplemental Figure 4).

In the whole series, a multivariate Cox proportional hazards model with FFS as the dependent variable and including the MRS and the IPS, we found that only the MRS was significant (Table 2). No interaction was observed between the IPS or the individual IPS variables and the MRS low- and high-risk groups (supplemental Table 5). Thus, to compare the molecular risk algorithm and the individual IPS components, a backward stepwise selection Cox model was tested, with FFS as the dependent variable and including the MRS and the individual components of the IPS in the global series of samples. Only MRS and stage IV were statistically significant (Table 2) and so were retained in the final Cox regression model. Patient stratification into quartiles on the basis of the Cox model identified a subgroup of advanced cHL patients (fourth quartile) with a very poor outcome: 5-year FFS of 24.3% (P < .001; Figure 2).

Here we describe, in a series of patients with advanced cHL, a 4-cluster/11-gene model derived from an initial selection of 30 potentially predictive markers that can be detected by RT-PCR and integrated into a molecular risk algorithm that can identify patient subgroups with very different probabilities of treatment failure. This approach is determined on the basis of reliable quantitative RT-PCR techniques applicable to paraffin-embedded diagnostic samples and follows similar approaches taken in breast cancer and other tumor types.^23,24  This benefits from previous expression profiling studies performed in cHL,¹⁰  and the improved knowledge about the role of the tumoral cell and the microenvironment in the pathogenesis and outcome of this disease.^13,25,26  This MRS, calculated in an initial estimation set, was confirmed in an independent set. Genes included in the score were selected on the basis of previous gene expression profiling data generated in independent sets of patients^11,27  and can be classified in 4 functional pathways.

The final 4-cluster/11-gene model can additionally incorporate one of the well-established clinical variables (stage IV), thus integrating the main molecular characteristics of the tumors related with treatment response and tumor burden estimation in a single scoring system. The multivariate Cox model indicates that most patients with stage IV cHL and with a high MRS (≥ 0.3) will have a very poor outcome, with 5-year FFS probability of 24.3% and OS probability of 76.3%. Therefore, this combination of stage IV and high MRS identifies a group of patients with very bad outcome who could have been initially missed by consideration of the IPS alone.^1,3 

It is of note that in the present series the IPS did not show any significant prognostic influence on FFS. From the individual components of the IPS, only stage IV disease remained significant in multivariate analyses. This finding is in agreement with recent studies in which the authors demonstrated that IPS is of limited utility in advanced HL cases treated in the modern era, where more accurate pathologic diagnosis, improved control of therapy, use of growth factors, and enhanced supportive care are yielding better outcomes compared with historic results.³  Thus, deviations from the standard therapy cannot be justified on the basis of the survival prediction capacity of the IPS, and identification of high-risk populations needs to be supplemented with molecular markers.

In the model, the expression of BCL2, BCL2L1, CASP3, HMMR, CENPF, CCNA2, CCNE2, and CDC2 included in the apoptosis and cell-cycle pathways, respectively, were correlated with short FFS. The expression of various antiapoptotic BCL2 family proteins has been repeatedly reported in HRS cells,^7,28,29  thereby contributing to the survival of HRS cells. BCL2 and BCL2L1 (BCL-XL) both frequently are expressed by HRS cells in cHL, and their levels have been associated with inferior FFS in patients treated with ABVD or equivalent regimens,⁵  thus confirming the importance of this group of apoptotic regulators for cHL outcome prediction. The second signature, cell cycle, is mainly composed of genes coding for regulatory proteins of the S and G₂/M phases of the cell cycle, thus directly related with cell proliferation. Again, expression of some of these markers has been previously described at the protein level in cHL,^7,30,31  and a significant prognostic value for CDC2 (CDK1) and CCNA2 (cyclin A2) protein expression in non-Hodgkin lymphomas was also found for both disease-free survival and OS. Interestingly, this aberrant association between increased expressions of antiapoptotic proteins and growth fraction-associated proteins in HRS cells provides further evidence that cell cycle and apoptosis regulation are profoundly disturbed and closely related in the disease, further justifying the inclusion of these 2 pathways in the predictive model.

Moreover, pathways involved in cell cycle and apoptosis regulation are rational therapeutic targets. Indeed, inhibitors of Cyclin-Cdk complexes (including Cyclin E-Cdk2, Cyclin A-Cdk2, and Cyclin B-Cdk1 complexes) are currently under preclinical and clinical investigation in different cancer types,^32-34  and these drugs could be considered for the treatment of advanced and refractory cHL patients. Likewise, the potential for targeting BCL2-related proteins in lymphoma is promising. Small molecule inhibitors of the Bcl-2 family have demonstrated high target affinity and an improved toxicity profile, and clinical trials of these agents are yielding interesting results.³⁵ 

Confirming previous observations about the importance of the reactive microenvironment for cHL patient outcome,^9,11,12 LYZ and STAT1 genes, expressed at high levels in a subset of tissue monocytes and activated macrophages, also are included in this model, and correlated with prolonged FFS and better outcome. The relevance of the cell composition of the reactive background in cHL has been reinforced by the data recently reported by Steidl et al.³⁶  They used gene expression profiling to identify a gene signature of tumor-associated macrophages that is associated with treatment failure, in an approach methodologically similar to previous reports from our group and others.^10,37  The discrepancy in the results concerning the role of macrophages may have arisen from technical differences (RT-PCR vs microarray gene-expression and immunohistochemistry) or the selection of markers such as LYZ and STAT1 in this study, reflecting a specific functional status of the monocyte-macrophages.³⁸ 

In addition, IRF4 (MUM1) expression was associated with longer FFS. This gene is an interferon regulatory factor, lymphocyte specific, induced after nuclear factor-κB (NF-κB) activation, which controls B-cell proliferation and differentiation and has recently been shown to be up-modulated by CD40 engagement in HL cells.³⁹  Interestingly, the lack of IRF4 protein has been previously associated with outcome in cHL,⁴⁰  representing a potential adverse prognostic factor. In addition, both IRF4 and BCL2L1 represent well-known NF-κB target genes^41,42  whose expression is induced after NF-κB pathway activation. Thus, our final model includes important subrogates from the NF-κB activation, which is thought to be an essential pathogenetic mechanism in this disease. It is of note that Bednarski et al⁴³  recently described how the inhibition of the canonical NF-κB pathway enhances the proapoptotic effects of adriamycin, thus also identifying NF-κB inhibition as an interesting therapeutic approach.

In conclusion, we have developed a molecular risk algorithm, on the basis of feasible and reproducible molecular techniques, that is capable of stratifying at the moment of diagnosis advanced cHL patients with different outcome. Moreover, a combination of this algorithm with the presence of clinical stage IV disease can be used to identify a group of cHL patients with a very bad outcome who could benefit from more intensive therapeutic approaches.

These results are promising, but further validation in larger and independent series of patient is needed for the model to become established as part of the clinical routine. Also the predictive value of this model should also be tested in patients treated with modern intensive chemotherapy.

We thank Nuria Malats, from the Genetic and Molecular Epidemiology Group at the CNIO, for her useful advice and supervision of the statistical analysis. We also thank Marién Castillo and Laura Cereceda, at the CNIO Tumor Bank, for collecting the human tumor samples and for their excellent assistance with data management. Members of the Spanish Hodgkin Lymphoma Study Group are cited in the supplemental Appendix.

This work was supported by grants from the Fondo de Investigaciones Sanitarias (PI08/1985, PI05/1623, PI05/2800, PI05/2327, RETIC RD06/0020/0107) and the Ministerio de Ciencia y Tecnología (SAF2008-03 871), Spain. B.S.-E. is supported by a grant from the Ministerio de Ciencia e Innovación (FIS), Spain.

Contribution: B.S.E., C.M., M.A.P., and J.F.G. contributed to the conception and design of the study, the analysis and interpretation of the data, and the drafting of the article; B.S.E. and A.L. performed statistical analysis; M.M.M. managed tissue banking; J.M., P.S., C.R.M., A.L., R.R., J.R., A.C., C.C., M.C., J.A., R.A., A.A., A.S., S.S., A.B., J.M.M., P.S.-G., F.B., C.R., M.F.F., J.G.L., M.G.-C., C.S., J.L.L., M.L.I., M.M., J.G.C., A.M., J.F., R.G.C, and J.F.T. managed patient databases and contributed with tumor samples and clinical follow-up; and all the authors read and approved the final version of the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Members of the Spanish Hodgkin Lymphoma Study Group are listed in the supplemental Appendix.

Correspondence: Miguel A Piris, MD, Lymphoma Group, Spanish National Cancer Centre (CNIO). E-28029 Madrid, Spain; e-mail: mapiris@cnio.es.