Abstract
Purpose: Plasma cell leukemia (PCL) is a rare form of plasma cell dyscrasia that presents either as a progression of previously diagnosed multiple myeloma, namely secondary PCL, or as initial manifestation of disease, namely primary PCL (pPCL). Although the presenting signs and symptoms include those seen in multiple myeloma, pPCL is characterized by several aspects that define a more aggressive course. Here, we have investigated the transcriptome of pPCLs and correlated differential expression profiles with outcome to provide insights into the biology of the disease.
Experimental Design: The expression profiles of 21 newly diagnosed pPCLs included in a multicenter prospective clinical trial were generated using high-density microarray, then evaluated in comparison with a representative series of patients with multiple myeloma and in association with clinical outcome.
Results: All but one of the pPCLs had one of the main immunoglobulin heavy-chain locus translocations, whose associated transcriptional signatures resembled those observed in multiple myeloma. A 503-gene signature distinguished pPCL from multiple myeloma, from which emerged 26 genes whose expression trend was associated with progressive stages of plasma cells dyscrasia in a large dataset from multiple institutions, including samples from normal donors throughout PCL. Finally, 3 genes were identified as having expression levels that correlated with response to the first-line treatment with lenalidomide/dexamethasone, whereas a 27-gene signature was associated with overall survival independently of molecular alterations, hematologic parameters, and renal function.
Conclusions: Overall, our data contribute to a fine dissection of pPCL and may provide novel insights into the molecular definition of patients with poorer prognosis. Clin Cancer Res; 19(12); 3247–58. ©2013 AACR.
pPCL is a very aggressive and rare hematologic malignancy. So far, genomic and clinical differences between pPCL and multiple myeloma have been shown, mainly on the basis of retrospective studies. Here, we took advantage of a prospective series of pPCLs included in a phase II clinical trial to investigate pPCL global mRNA expression profiles. A 503-gene signature that distinguished pPCL from multiple myeloma was identified; a fraction of these genes have been suggested whose expression is putatively associated with the aggressiveness of clinical presentation. In addition, a 27-gene model has been identified with potential clinical relevance that distinguished, within pPCL, those with worst outcome. Because pPCL represents a high-risk clinical entity per se, our findings help provide insights into molecular features that could prognostically stratify patients with this aggressive form of plasma cell dyscrasia.
Introduction
Plasma cell leukemia (PCL) is a rare but highly aggressive disease that represents 3% to 5% of all plasma cell (PC) disorders. Its diagnosis is based on the Kyle criteria, which requires circulating PCs to account for at least 20% of peripheral blood leukocytes and/or an absolute circulating PC count of 2 × 109/L, with evidence of monoclonal gammopathy. Clinically, the PCL can be distinguished into primary PCL (pPCL), originating de novo without any prior history of multiple myeloma, or secondary PCL (sPCL), arising from a preexisting multiple myeloma tumor that eventually progressed to the leukemic phase. The prognosis of pPCL is very poor and sPCL has the worst prognosis, with median survival durations of 7 to 11 months and 2 to 7 months, respectively. The pathogenetic mechanisms underlying either the de novo presentation or the aggressive transformation of relapse/refractory multiple myeloma are largely unknown (1).
So far, few studies on both pPCL and sPCL have been conducted with an attempt to characterize the molecular alterations underlying PCL biology, mainly on the basis of immunophenotypic, fluorescence in situ hybridization (FISH) and cytogenetic or array-based comparative genomic hybridization (aCGH) analyses (2–6). All of them, albeit suffering restrictions of numerically limited cohorts of cases, agreed in identifying significant molecular differences among the multiple myeloma and PCL presentation, specifically in regard to the surface immunophenotype (the differential expression of CD20, CD56, CD9, CD117, and HLA-DR antigens) and the increased number of genomic aberrations identified in PCL (7). Indeed, PCL shows higher incidence of translocations involving immunoglobulin heavy-chain locus (IGH@) at 14q32, in particular t(11;14), whereas very few patients with hyperdiploidy have been reported differently from multiple myeloma (1). Among PCL, only t(4;14) was found to be associated with shorter overall survival (3), while controversial data exist on the favorable prognostic significance of t(11;14) [1]. MYC rearrangement, a late genetic event correlated with disease progression in multiple myeloma, has been frequently found in PCL also in association with a worse trend of the disease (5, 6).
To date, high-resolution transcriptional profiling studies of PCL are limited to a retrospective small series reported by us (8) and the recently published investigation from Usmani and colleagues (9), who described pPCL cases included in the Total Therapy protocols database. Both studies indicated that transcriptional features may distinguish the 2 entities; moreover, in the latter study the authors additionally indicated that pPCL is an independent adverse prognostic factor in comparison with multiple myeloma.
Here, we investigated a series of 21 previously untreated pPCL patients included in a prospective clinical trial using whole-transcript profiling arrays, and correlated the transcriptional profiles with the primary and secondary outcome endpoints to evaluate whether different molecular or clinical entities may exist within pPCL. In addition, we compared the transcriptional profiles of the 21 pPCL with those of a representative dataset of 55 newly diagnosed primary multiple myeloma tumors; the differentially expressed genes were then evaluated in a large dataset, profiled on a different generation array platform, from multiple institutions including normal donors, different forms of PC dyscrasia from asymptomatic to aggressive stages, and a collection of human myeloma cell lines, with the aim to identify transcripts that might contribute to the pathogenesis of pPCL as well as to tumor progression in PC dyscrasia.
Materials and Methods
Patients
Pathologic specimens were collected from 23 untreated pPCL patients included in a multicenter Italian clinical trial (RV-PCL-PI-350, EudraCT N° 2008-003246-28), an open-label, exploratory, single-arm, 2-stage study aimed at evaluating safety and antitumor activity of combined lenalidomide and dexamethasone (LD) as the first-line treatment of previously untreated pPCLs. All patients gave their informed consent for molecular analyses. The primary endpoint of the study was the response rate after a 4-cycle therapy over a 4-month schedule, according to International Uniform Criteria (10, 11); secondary endpoints were: (i) time to progression; (ii) progression-free and overall survival (OS); (iii) eligibility to undergo autologous or allogeneic stem cells transplantation after LD treatment; and (iv) serious/severe adverse event rate. Until February 2012, the median follow-up was 23 months (range, 9–32). To the purpose of the present work, we considered response rate and OS. Highly purified (≥90%) bone marrow PCs were obtained as previously reported (8). Demographic information and details on molecular alterations of the whole cohort were reported in a previous manuscript aimed at the genomic characterization of the samples (12). Fifty-five newly diagnosed multiple myeloma patients were included in the study; these cases have been previously reported (13) and were selected on the basis of the representativeness of the main molecular characteristics.
Generation of gene expression data
Total RNA was available from 21 pPCL and 55 multiple myeloma patients. Multiple myeloma samples were part of a larger cohorts previously profiled by us on old-generation array (GeneChip HG-U133A; ref. 13), herein resampled on newer generation arrays and selected on the basis of their representativeness of molecular features, to allow comparison with pPCL profiles. Preparation of DNA single-stranded sense target, hybridization to GeneChip Gene 1.0 ST arrays (Affymetrix Inc.), and scanning of the arrays (7G Scanner, Affymetrix Inc.) were carried out according to the manufacturer's instructions. Log2-transformed expression values were extracted from CEL files and normalized using NetAffx Transcript Cluster Annotations, Release 32 (June 2011) and Robust Multi-array Average (RMA) procedure in Expression Console software (Affymetrix Inc.). Non-annotated transcript clusters were discarded. The expression values of transcript cluster ID specific for loci representing naturally occurring read-through transcriptions or mapped to more than 1 chromosomal location were summarized as median value for each sample. The data have been deposited at NCBI Gene Expression Omnibus repository (http://www.ncbi.nlm.mih.gov/geo) and are accessible through GEO Series accession number GSE39925. The multiple myeloma samples were stratified into the 5 translocation/cyclin groups as previously described (14).
The dataset of samples profiled on GeneChip HG-U133A arrays (Affymetrix Inc.) was generated by 2 proprietary (GSE13591 and GSE6205) and 2 publicly available datasets (GSE6477 and GSE6691; refs. 15, 16), which globally included 24 healthy donors (N), 33 monoclonal gammopathy of undetermined significance (MGUS), 24 smoldering multiple myeloma (SMM), 200 newly diagnosed and 26 relapsed multiple myeloma, 9 PCL patients (3 of whom primary), together with 23 HMCLs fully characterized in a previous study (17). Expression values were extracted from CEL intensity files using GeneAnnot custom chip definition files version 2.2.0 and normalized by Robust Multi-array Average procedure, as previously described (18). To prevent batch effects, namely the inclusion of low-quality or nonreproducible data, normalized unscaled standard error (NUSE) and relative log-expression (RLE) distributions were generated in aroma.affymetrix package for all samples; samples were removed if the 25th or the 75th percentile of NUSE and RLE exceeded the value of ±1.05 or ±0.5, respectively, which led to a definitive amount of 323 comparable samples.
Microarray data analysis
Hierarchical agglomerative clustering was carried out on the list of genes whose average change in expression levels varied at least 2-fold from the mean across the considered data set, as previously reported (8). Euclidean and Ward were respectively used as distance and linkage methods in hclust function in R software.
The differentially expressed genes discriminating the considered classes were identified using the Significant Analysis of Microarrays software version 4.00 as previously described (19). The functional annotation analysis on the selected lists was carried out using the gene set enrichment analysis (GSEA; ref. 20) to identify significantly enriched curated gene sets (version 3.0) among those included in the MSigDB database (http://www.broadinstitute.org/gsea/msigdb/index.jsp). A false discovery rate that was lower than 25% was used as an enrichment statistic, and the rank was permuted with 1,000 permutations. Only gene sets including 5 through 1,000 genes were considered in the analysis.
The Jonckheere–Terpstra test was used with the clinfun package in R software to investigate the significance of the trend from normal donors through PCL cases in the expression levels of the selected gene list generated on HG-U133A arrays. To avoid biases occurring conceivably due to variability in the target detection technology, the genes whose trend in multiple myeloma and PCL cases on HG-U133A arrays was incoherent with the data generated on HuGene1.0st arrays were not considered. For robustness, to reduce biases that might be due to numerical imbalances within the groups (with the multiple myeloma group being largely overrepresented) and to gain the most significant trends, an additional criterion was imposed that only steadily ascending or descending median values were allowed; in addition, a differential expression should be observed in at least 1 condition (Kruskal–Wallis tests conducted using the stat package in R software). The Benjamini and Hochberg correction was used to adjust significance of multiple tests.
Survival analysis
A Cox proportional hazards model and 100.000 permutations were used in the globaltest function of R software to search for a correlation between the gene expression levels and OS (21). Patients were stratified into 2 groups using the prognostic gene signatures and hierarchical clustering with Euclidean and Ward as distance and linkage metrics, respectively. The groups identified by this approach were tested for association with survival using the Kaplan–Meier estimator and log-rank test, and P-values were calculated according to the standard normal asymptotic distribution. Survival analysis was conducted with the survcomp package in R software. Independence between commonly used prognostic factors and gene signatures was tested using multivariate Cox regression procedure of survival R package.
Real-time quantitative polymerase chain reaction
The expression levels of NFKBIA, RELB, RHOA, and SDC1 were analyzed in purified CD138+ cells by means of real-time quantitative PCR (qRT-PCR). Total RNA was converted to cDNA using M-MLV reverse transcriptase (Invitrogen). Inventoried or Made-to-Order TaqMan gene expression assays (Hs00153283_m1 for NFKBIA, Hs00232399_m1 for RELB, Hs01051295_m1 for RHOA, and Hs00896423_m1 for SDC1) and the TaqMan Fast Universal PCR Master Mix were used according to the manufacturer's instructions (Applied Biosystems). GAPDH-specific predeveloped assay reagent (PDAR, Applied Biosystems) was used as the internal control. The measurement of gene expression was done using the ABI PRISM 7900HT Sequence Detection System (Applied Biosystems). All the samples were run in duplicate. Data were expressed as 2−ΔCt (Applied User Bulletin No. 2).
Results
Association between primary IGH translocation and the transcriptional profiles of pPCL tumors
Here, we analyzed the transcriptional profile of 21 of the 23 pPCLs (for whom RNA material was available) included in the prospective Italian multicentric clinical trial RV-PCL-PI-350, all of which were characterized by FISH for the main genomic aberrations. Specifically, 3 of them had t(4;14), 9 had t(11;14), 7 had t(14;16), and 1 had t(14;20); notably, only 1 patient lacked any of these major IGH translocations (12). In addition, we generated the gene expression profiles of 55 primary multiple myeloma tumors, representative of the type and frequency of the major molecular aberrations of this clinical entity. Unsupervised analysis was carried out on the 1,145 most variable transcripts in the 21 pPCLs dataset to evaluate whether natural grouping of the transcriptional profiles might be associated with distinct molecular subgroups. Indeed, as shown in Fig. 1A, the hierarchical clustering was mostly driven by the presence of the main recurrent IGH chromosomal translocations. Specifically, all cases with translocated MAF genes clustered in 1 main branch (P = 4.423e-05), whereas the t(11;14)- and t(4;14)-positive pPCLs were clustered together in the second branch. Within the MAF-translocated PCL, a sample was clustered lacking major translocations (PCL-020), with expression profile resembling those of MAF/MAFB samples, in all likelihood due to the exclusive and spiked expression of MAFA, another member of the MAF family transcription factor (Supplementary Fig. S1). Unfortunately, lack of available material did not allow investigation of the occurrence of a translocation event.
A, hierarchical clustering analysis of gene expression profiles of 21 pPCL cases. Samples are grouped according to the expression levels of the 1,145 most variable genes. Main molecular alterations are shown; black indicates occurrence. B, heatmap of the 47 differentially expressed genes identified by multi-class analysis of 21 pPCL patients stratified into the 3 main IGH translocation groups. Visualization of the expression in multiple myeloma samples is included. In the heatmap, the color-scale bar represents the relative gene expression changes normalized by the SD, and the color changes in each row represent gene expression relative to the mean across the samples.
A, hierarchical clustering analysis of gene expression profiles of 21 pPCL cases. Samples are grouped according to the expression levels of the 1,145 most variable genes. Main molecular alterations are shown; black indicates occurrence. B, heatmap of the 47 differentially expressed genes identified by multi-class analysis of 21 pPCL patients stratified into the 3 main IGH translocation groups. Visualization of the expression in multiple myeloma samples is included. In the heatmap, the color-scale bar represents the relative gene expression changes normalized by the SD, and the color changes in each row represent gene expression relative to the mean across the samples.
We conducted a multi-class supervised analysis to identify genes that might specifically differentiate pPCL cases stratified on the main IGH chromosomal translocations. Forty-seven genes were differentially expressed in the 3 classes (Fig. 1B; Supplementary Table S1). As expected, the direct targets of t(11;14) (CCND1), t(4;14) (WHSC1), and MAF genes were found to be respectively upregulated in each corresponding pPCL subgroup (to the aim of the analysis, PCL-020 was considered part of the MAF group). The FGFR3 gene was not identified in the gene signature of the t(4;14), likely due to a missing derivative chromosome 14 in 1 of the 3 cases. Notably, the CCND2 transcript showed the highest score and positively marked both the MAF cases and, although at a lesser extent, the t(4;14)-translocated pPCLs. As depicted in Fig. 1B, the signatures of the IGH translocations largely resembled those observed in multiple myeloma and, in part, previously described.
Identification of a 503-gene transcriptional signature distinguishing pPCL from multiple myeloma
On the basis of its clinical presentation, it is conceivable that pPCL is characterized by distinct transcriptional patterns from intramedullary multiple myeloma. We, therefore, investigated the relationship between the transcriptional profiles of pPCL and multiple myeloma tumors. First, we carried out hierarchical clustering of the 1,166 most variable genes in the dataset, including the 21 pPCL and 55 multiple myeloma cases (Fig. 2). Again, the clustering was mainly driven by the IGH translocations: in fact, t(14;16)/t(14;20) or t(4;14) or t(11;14) cases clustered together independent of belonging to multiple myeloma or pPCL groups, suggesting that such translocation events had stronger consequences on the transcriptional fingerprint than those due to leukemic phenotype.
Hierarchical clustering analysis of gene expression profiles of 21 pPCL and 55 multiple myeloma cases. Samples are grouped according to the expression levels of the 1,166 most variable genes. Main molecular alterations are indicated as black boxes, as well as the pPCL type in the last lane.
Hierarchical clustering analysis of gene expression profiles of 21 pPCL and 55 multiple myeloma cases. Samples are grouped according to the expression levels of the 1,166 most variable genes. Main molecular alterations are indicated as black boxes, as well as the pPCL type in the last lane.
Next, in order to identify the genes that specifically distinguished the 2 clinical entities, we conducted a supervised analysis between the 21 pPCL and the 55 multiple myeloma cases. We found 503 differentially expressed genes, 366 of which upregulated in pPCLs (Supplementary Table S2). Functional annotation analysis was conducted using GSEA (Supplementary Table S3). Among the identified gene sets, it is worth mentioning those associated with the NF-kB pathway, structural organization of the cell, and migration, together with some gene sets that suggested putative alterations of CD40, TGFB, AKT, and FAS signaling pathways. Notably, the SDC1 gene, codifying the CD138 surface antigen, showed significantly lower expression level in pPCLs. Prompted by these findings, on the basis of gene expression data, we evaluated whether pPCL cases were likely to show more activated NF-kB pathway or higher proliferation pattern. To this aim, we took advantage of the NF-kB indexes on the basis of expression data as described by Annunziata and colleagues (22) and Keats and colleagues (23); and the proliferation indexes as described by Hose and colleagues (24), Bergsagel and colleagues (25), and Zhan and colleagues (26). Interestingly, both the procedures for NF-kB indexes calculation highlighted increased NF-kB pathway activity in pPCL (Supplementary Fig. S2A and 2B); conversely, no significant difference was raised according to the gene proliferation indexes (Supplementary Fig. S2C–2E). As a methodologic validation of our results, 4 of the genes identified and included in the mentioned gene sets (NFKBIA, RELB, RHOA, and SDC1) were validated by qRT-PCR in all samples for which RNA was available (15 pPCL and 40 multiple myeloma patients; Supplementary Fig. S3). The correlation coefficients of the expression of each gene as determined by microarray or qRT-PCR indicated an optimal concordance for all genes (median Pearson correlation coefficient, 0.73).
Modulation of genes included in multiple myeloma-pPCL signature throughout the different stages of PC dyscrasia
We hypothesized that the identified 503-gene signature distinguishing pPCL and multiple myeloma cases could be accountable to the occurrence of overt/aggressive disease. Therefore, we tried to elucidate whether the genes of this list, or at least a part of them, may be modulated through (if not discriminate) the different stages of the disease; namely, whether the expression trend of these transcripts could underlie the transition from indolent forms to aggressive stages. To verify this, we generated a comprehensive dataset including independent samples from 3 different cohorts, encompassing PC samples from normal subjects, MGUS, smoldering multiple myeloma, newly diagnosed multiple myeloma, and PCL patients, all profiled on the HG-U133A microarray chip. Of the public available samples, 274 passed the quality controls applied to prevent biases due to inter-cohort analysis. In this dataset, we investigated the expression profiles of 360 of the 503 differentially expressed transcripts represented on the HG-U133A array, and tested whether a trend existed in their expression levels through all the progressively malignant stages of PC dyscrasia. We selected the 26 most significant genes (Supplementary Fig. S4A) whose trend in expression levels correlated with progressive disease. Larger expression level spreads were observed in multiple myeloma subgroups than in other stages, in line with the well-accepted notion that multiple myeloma represents a widely heterogeneous entity. However, no significant associations were found between the expression levels of the 26 genes and the molecular stratification of the myeloma samples on the basis of the occurrence of IGH translocations, with few marginally significant exceptions (Supplementary Fig. S4B). In addition, we evaluated the expression profiles of those 360 genes in 26 relapsed multiple myeloma and 23 HMCLs. Notably, although not included in the analysis, the trend of the expression values of all the genes in HMCLs was concordant with that observed in primary tumors, thus reinforcing the suggestion that the deregulation of these genes (either with increasing or decreasing pattern) is compatible with the “strength” of neoplastic transformation. Among the identified transcripts, it is worth mentioning the proteasome-associated gene PSMD6 and the polycomb group protein codifying gene EZH2, whose expression has been associated with tumor burden in multiple myeloma (27).
Identification of transcriptional profiles with clinical relevance in the definition of pPCL with poorer outcome
Finally, the transcriptional features of pPCL were evaluated in the context of the outcome data of the prospective series of pPCL patients. In particular, we evaluated whether the occurrence of specific gene signature might be associated either with response rate or OS. To this aim, we first investigated the dataset for differentially expressed genes in patients who failed to respond to frontline 4-cycle therapy of LD. The analysis led to the identification of 3 genes (YIPF6, EDEM3, and CYB5D2) able to distinguish nonresponder patients from responders [i.e., complete response (CR), very good partial response (VGPR), and partial response (PR); Fig. 3]. When PR, VGPR, and CR were considered separately, no specific differentially expressed genes could be evidenced. The identified transcriptional pattern will be helpful to integrate results on efficacy and side effects of the first-line treatment of LD in pPCL.
Boxplot distribution of the expression levels of the 3 differentially expressed genes identified by supervised analysis of pPCL patients stratified according to response rate. NR, nonresponder; PR, partial remission; VGPR, very good partial remission; CR, complete remission.
Boxplot distribution of the expression levels of the 3 differentially expressed genes identified by supervised analysis of pPCL patients stratified according to response rate. NR, nonresponder; PR, partial remission; VGPR, very good partial remission; CR, complete remission.
Next, we assessed the relationship between each of the 1145 most variable genes across the pPCL dataset and OS, using a statistical approach based on the linear regression model in which the distribution of the response variable is modeled as a function of the expression levels of each gene. Of the 1145 genes, 27 reached a highly significant correlation (P < 0.01) with OS (Table 1). On the basis of this 27-gene model, the 18 pPCL cases for whom follow-up information was available could be divided into 2 groups, of 6 and 12 patients respectively, who showed different outcomes (Fig. 4A and B). This model retained independency from all the molecular characteristics available (Table 2), as well as from age, sex, lactate dehydrogenase (LDH) levels, renal function, and hematologic parameters (data not shown). However, it is worth mentioning that none of the cytogenetic aberration was associated per se with OS (12). It is also worth emphasizing that the 27-gene model was not independent of patients being subjected to autologous stem cell transplantation (ASCT), indicating that this therapeutic approach points definitively toward a more favorable outcome.
A, heatmap of the pPCL samples clustered on the 27 genes identified by globaltest as significantly associated with OS. Molecular features are indicated above the heatmap (black indicates occurrence). B, Kaplan–Meier curves of the 2 groups defined by the 27-gene model. The colors of the groups correspond to those shown in the dendrogram in A.
A, heatmap of the pPCL samples clustered on the 27 genes identified by globaltest as significantly associated with OS. Molecular features are indicated above the heatmap (black indicates occurrence). B, Kaplan–Meier curves of the 2 groups defined by the 27-gene model. The colors of the groups correspond to those shown in the dendrogram in A.
List of the 27 genes with P < 0.01 in globaltest analysis of expression data and OS
Gene symbol . | Gene title . | Cytoband . | P-value . | Corr with survival . |
---|---|---|---|---|
PECAM1 | Platelet/endothelial cell adhesion molecule | 17q23.3 | 0.0006 | pos |
MKX | Mohawk homeobox | 10p12.1 | 0.0020 | pos |
FAM111B | Family with sequence similarity 111, member B | 11q12.1 | 0.0021 | neg |
MCTP1 | Multiple C2 domains, transmembrane 1 | 5q15 | 0.0022 | neg |
CALCRL | Calcitonin receptor-like | 2q32.1 | 0.0023 | pos |
C10orf10 | Chromosome 10 open reading frame 10 | 10q11.21 | 0.0028 | neg |
FNBP1 | Formin binding protein 1 | 9q34.11 | 0.0028 | neg |
EFEMP1 | EGF-containing fibulin-like extracellular matrix protein 1 | 2p16.1 | 0.0030 | neg |
C3orf14 | Chromosome 3 open reading frame 14 | 3p14.2 | 0.0031 | pos |
ALDH1L2 | Aldehyde dehydrogenase 1 family, member L2 | 12q23.3 | 0.0032 | pos |
WARS | Tryptophanyl-tRNA synthetase | 14q32.2 | 0.0034 | pos |
SLC15A2 | Solute carrier family 15 (H+/peptide transporter), member 2 | 3q13.33 | 0.0043 | pos |
FAIM3 | Fas apoptotic inhibitory molecule 3 | 1q32.1 | 0.0043 | neg |
CPEB4 | Cytoplasmic polyadenylation element binding protein 4 | 5q35.2 | 0.0048 | neg |
EDN1 | Endothelin 1 | 6p24.1 | 0.0055 | neg |
PVALB | Parvalbumin | 22q12.3 | 0.0056 | neg |
LY86 | Lymphocyte antigen 86 | 6p25.1 | 0.0058 | neg |
LAPTM5 | Lysosomal protein transmembrane 5 | 1p35.2 | 0.0065 | neg |
RNU5D | RNA, U5D small nuclear | 1p34.1 | 0.0077 | pos |
PARP15 | Poly (ADP-ribose) polymerase family, member 15 | 3q21.1 | 0.0079 | neg |
PLEKHF2 | Pleckstrin homology domain containing, family F (with FYVE domain) member 2 | 8q22.1 | 0.0085 | neg |
PDK4 | Pyruvate dehydrogenase kinase, isozyme 4 | 7q21.3 | 0.0089 | neg |
TNFAIP3 | Tumor necrosis factor, alpha-induced protein 3 | 6q23.3 | 0.0089 | neg |
FAM105A | Family with sequence similarity 105, member A | 5p15.2 | 0.0091 | neg |
CTH | Cystathionase (cystathionine gamma-lyase) | 1p31.1 | 0.0091 | pos |
HOOK1 | Hook homolog 1 (Drosophila) | 1p32.1 | 0.0092 | pos |
TCN2 | Transcobalamin II | 22q12.2 | 0.0094 | neg |
Gene symbol . | Gene title . | Cytoband . | P-value . | Corr with survival . |
---|---|---|---|---|
PECAM1 | Platelet/endothelial cell adhesion molecule | 17q23.3 | 0.0006 | pos |
MKX | Mohawk homeobox | 10p12.1 | 0.0020 | pos |
FAM111B | Family with sequence similarity 111, member B | 11q12.1 | 0.0021 | neg |
MCTP1 | Multiple C2 domains, transmembrane 1 | 5q15 | 0.0022 | neg |
CALCRL | Calcitonin receptor-like | 2q32.1 | 0.0023 | pos |
C10orf10 | Chromosome 10 open reading frame 10 | 10q11.21 | 0.0028 | neg |
FNBP1 | Formin binding protein 1 | 9q34.11 | 0.0028 | neg |
EFEMP1 | EGF-containing fibulin-like extracellular matrix protein 1 | 2p16.1 | 0.0030 | neg |
C3orf14 | Chromosome 3 open reading frame 14 | 3p14.2 | 0.0031 | pos |
ALDH1L2 | Aldehyde dehydrogenase 1 family, member L2 | 12q23.3 | 0.0032 | pos |
WARS | Tryptophanyl-tRNA synthetase | 14q32.2 | 0.0034 | pos |
SLC15A2 | Solute carrier family 15 (H+/peptide transporter), member 2 | 3q13.33 | 0.0043 | pos |
FAIM3 | Fas apoptotic inhibitory molecule 3 | 1q32.1 | 0.0043 | neg |
CPEB4 | Cytoplasmic polyadenylation element binding protein 4 | 5q35.2 | 0.0048 | neg |
EDN1 | Endothelin 1 | 6p24.1 | 0.0055 | neg |
PVALB | Parvalbumin | 22q12.3 | 0.0056 | neg |
LY86 | Lymphocyte antigen 86 | 6p25.1 | 0.0058 | neg |
LAPTM5 | Lysosomal protein transmembrane 5 | 1p35.2 | 0.0065 | neg |
RNU5D | RNA, U5D small nuclear | 1p34.1 | 0.0077 | pos |
PARP15 | Poly (ADP-ribose) polymerase family, member 15 | 3q21.1 | 0.0079 | neg |
PLEKHF2 | Pleckstrin homology domain containing, family F (with FYVE domain) member 2 | 8q22.1 | 0.0085 | neg |
PDK4 | Pyruvate dehydrogenase kinase, isozyme 4 | 7q21.3 | 0.0089 | neg |
TNFAIP3 | Tumor necrosis factor, alpha-induced protein 3 | 6q23.3 | 0.0089 | neg |
FAM105A | Family with sequence similarity 105, member A | 5p15.2 | 0.0091 | neg |
CTH | Cystathionase (cystathionine gamma-lyase) | 1p31.1 | 0.0091 | pos |
HOOK1 | Hook homolog 1 (Drosophila) | 1p32.1 | 0.0092 | pos |
TCN2 | Transcobalamin II | 22q12.2 | 0.0094 | neg |
Multivariate analysis comparing the 27-gene model with molecular variables and UAMS, IFM, and UK gene-risk models in pPCL series
Covariates . | HRa . | lo95%CIb . | up95%CIb . | P . |
---|---|---|---|---|
27-gene model | 33.28 | 2.77 | 400.09 | 5.73E-03 |
del(13q) | 4.89 | 0.37 | 64.92 | 2.29E-01 |
27-gene model | 20.54 | 2.42 | 174.38 | 5.61E-03 |
del(17p) | 0.93 | 0.19 | 4.41 | 9.24E-01 |
27-gene model | 36.37 | 2.76 | 478.67 | 6.27E-03 |
gain(1q) | 0.37 | 0.05 | 2.84 | 3.39E-01 |
27-gene model | 19.07 | 2.23 | 163.11 | 7.10E-03 |
del(1p) | 0.97 | 0.2 | 4.66 | 9.72E-01 |
27-gene model | 35 | 3.19 | 383.61 | 3.61E-03 |
t(11;14) | 3.2 | 0.47 | 22 | 2.36E-01 |
27-gene model | 20.26 | 2.35 | 174.38 | 6.15E-03 |
MAF-translocations | 1.28 | 0.26 | 6.41 | 7.65E-01 |
27-gene model | 16.28 | 1.79 | 148.15 | 1.33E-02 |
del(8p) | 0.83 | 0.14 | 4.81 | 8.33E-01 |
27-gene model | 16.30 | 1.68 | 158.46 | 1.61E-02 |
UAMS 70-gene model | 1.87 | 0.34 | 10.21 | 4.70E-01 |
27-gene model | 19.92 | 2.30 | 172.94 | 6.66E-03 |
UAMS 17-gene model | 1.24 | 0.22 | 6.92 | 8.08E-01 |
27-gene model | 20.54 | 2.43 | 173.92 | 5.56E-03 |
IFM 15-gene model | 0.79 | 0.15 | 4.20 | 7.79E-01 |
27-gene model | 20.63 | 2.40 | 177.53 | 5.84E-03 |
UK 6-gene model | 1.38 | 0.29 | 6.51 | 6.83E-01 |
Covariates . | HRa . | lo95%CIb . | up95%CIb . | P . |
---|---|---|---|---|
27-gene model | 33.28 | 2.77 | 400.09 | 5.73E-03 |
del(13q) | 4.89 | 0.37 | 64.92 | 2.29E-01 |
27-gene model | 20.54 | 2.42 | 174.38 | 5.61E-03 |
del(17p) | 0.93 | 0.19 | 4.41 | 9.24E-01 |
27-gene model | 36.37 | 2.76 | 478.67 | 6.27E-03 |
gain(1q) | 0.37 | 0.05 | 2.84 | 3.39E-01 |
27-gene model | 19.07 | 2.23 | 163.11 | 7.10E-03 |
del(1p) | 0.97 | 0.2 | 4.66 | 9.72E-01 |
27-gene model | 35 | 3.19 | 383.61 | 3.61E-03 |
t(11;14) | 3.2 | 0.47 | 22 | 2.36E-01 |
27-gene model | 20.26 | 2.35 | 174.38 | 6.15E-03 |
MAF-translocations | 1.28 | 0.26 | 6.41 | 7.65E-01 |
27-gene model | 16.28 | 1.79 | 148.15 | 1.33E-02 |
del(8p) | 0.83 | 0.14 | 4.81 | 8.33E-01 |
27-gene model | 16.30 | 1.68 | 158.46 | 1.61E-02 |
UAMS 70-gene model | 1.87 | 0.34 | 10.21 | 4.70E-01 |
27-gene model | 19.92 | 2.30 | 172.94 | 6.66E-03 |
UAMS 17-gene model | 1.24 | 0.22 | 6.92 | 8.08E-01 |
27-gene model | 20.54 | 2.43 | 173.92 | 5.56E-03 |
IFM 15-gene model | 0.79 | 0.15 | 4.20 | 7.79E-01 |
27-gene model | 20.63 | 2.40 | 177.53 | 5.84E-03 |
UK 6-gene model | 1.38 | 0.29 | 6.51 | 6.83E-01 |
NOTE: The model was available for all covariates with the exception of t(4;14), due to the low number of occurrence, and of ASCT, whose occurrence was inversely correlated with being part of the high-risk group that prevented the calculation of hazard ratio model in relationship with the configuration of survival data versus value of covariates.
aHR, hazard ratio; blo95%CI, up95%CI, lower and upper limits, respectively, of the 95% confidence interval of the hazard ratio.
We finally tested the independency of the 27-gene model from other gene-risk models on the basis of gene expression data described in multiple myeloma (28–30). Specifically, to this aim, we first stratified the samples included in our dataset according to the criteria defined by the UAMS 70- and 17-gene models from Shaughnessy and colleagues (30), the IFM 15-gene model by Decaux and colleagues (28), and the UK 6-gene model by Dickens and colleagues (29). Of them, only UAMS models effectively stratified pPCL patients into 2 significantly different risk groups (data not shown). Importantly, the 27-gene model retained significant correlation with outcome against both UAMS and other models (Table 2).
Discussion
In the present study, we reported a comprehensive molecular and transcriptional analysis of a series of pPCLs patients included in a prospective clinical trial aimed at exploring efficacy and safety of the LD combination as first-line therapy in previously untreated patients. The analysis was particularly aimed at identifying those transcriptional features that may contribute to explain the aggressive phenotype of PCL, as well as at correlating their gene expression profiles with the clinical outcome.
According to previous reports (2, 4–6), a high incidence of the 14q32 translocations was found in our series of pPCLs. The transcriptional analysis further supported that the main IGH chromosomal translocations drive pPCL patients clustering and are associated with specific signatures, as it occurs in multiple myeloma. The combined analysis of pPCL and multiple myeloma patients confirmed that pPCLs and multiple myelomas clustered together in specific IGH translocation groups, thus suggesting that the transcriptional effects related to these translocations have an impact that overcomes the signatures related to the specific disease. Moreover, the supervised analysis revealed that IGH translocations have transcriptional effects in pPCL cases that actually resembled, in a virtually identical fashion, those observed in multiple myeloma, identifying transcripts that have been shown as direct or indirect targets of the translocations in multiple myeloma, as previously described (14). Conversely, no association has been found in our dataset between cytogenetic abnormalities and clinical outcome. Therefore, occurrence of IGH translocations in a very large fraction, if not almost the totality, of pPCL cases was not conceivable to justify per se the more aggressive phenotype of leukemic forms than the large majority of multiple myeloma.
We thus aimed our analysis at investigating the transcriptional pathways differentially expressed in pPCL with regard to multiple myeloma condition. Of the 503-gene signature identified, several genes have been identified as related to cytoskeleton functions and Rho protein signaling pathway (31–33), involved in cell adhesion/migration processes. Our analysis also highlighted the involvement of NF-kB pathway-associated genes, specifically acting as functional regulators in NF-kB pathway or potentially responsive to NF-kB transcriptional modulation. In particular, the RelB subunit involved in NF-kB noncanonical pathway, whose activation was found to be associated to the cell adhesion-mediated drug resistance (CAM-DR) in multiple myeloma (34), was upregulated in pPCLs as well as the inhibitory subunits NFKBIE (IkBϵ) and NFKBIA (IkBα). In addition, TRAF1, TRAF4, and TRAF5, members of the family of TNF–receptor-associated factor with different roles in NF-kB pathway, were expressed at higher level in pPCLs (35, 36). Overall, these considerations, together with the concurrent findings of higher NF-kB indexes in pPCL cases, are suggestive and compatible with more activated NF-kB pathway and aggressive phenotype in pPCL. Conversely, a limited, if not absent, discrepancy between pPCL and multiple myeloma has been identified, according to published expression-based indexes (24–26), as regards proliferation. This might be explained in that such indexes were mainly derived by the comparison between primary tumors and cell lines that might, therefore, account for more enhanced proliferation activity in relapsed/refractory myeloma, evolving to a more aggressive and uncontrolled leukemic progression. The transcriptional pattern of pPCL was recently investigated by Usmani and colleagues, who characterized the clinicobiological features of a series of pPCL patients, partly included in the Total Therapy protocols trials (9), suggesting that pPCL represents a highly significant adverse feature of myeloma tumors in relationship with both OS and PFS as well as CR duration. This was further supported by the occurrence of high-risk variables (advanced age; abnormal albumin, LDH, and β2-microglobulin levels; chromosome abnormalities, etc.) that undoubtedly confer to pPCL a highly aggressive phenotype with very poor outcome. They described a 203-gene signature, obtained from 20 pPCL cases, that only in a minimal part (∼15%) overlapped with our 503-gene list. As possible explanations to this finding, the new-generation arrays used in our analysis (leading to only 125 transcripts commonly represented on the 2 platforms) and the different balance between pPCL and non-pPCL samples in the 2 datasets could be considered. Moreover, the authors did not provide specific information on the molecular alterations (translocations, deletions, and amplifications) of their cases; therefore, we could not exclude that different stratification within pPCL cases might affect the differential expression analysis. Furthermore, their data were not available in public repositories, which prevented any meta-analyses, validation, or verification of similarities and differences between the 2 datasets.
On the basis of the list of differentially expressed genes between pPCL and multiple myeloma, we took advantage of a composite dataset from multiple institutions to investigate whether the modulation of expression levels of some transcripts included in the 503-gene signature might be correlated with progression in PC dyscrasias, provided that the increased or decreased expression could be due to or associated with features characterizing highly injured PCs and/or aggressive clinical course. Among the 26 genes with a significant trend (8 decreasing and 18 increasing) from normal PC condition to PCL through the different forms of PC dyscrasia, it is worth mentioning the methyltransferase EZH2, whose expression correlates with tumor burden during multiple myeloma disease progression and which is constitutive in IL–6-independent cell lines (27), and that was among the most overexpressed genes associated with the proliferation (PR) multiple myeloma subgroup in the UAMS classification (26). Overall, the 26 genes identified here showed significant modulation in the expression levels through the different progressive presentation of PC dyscrasia and exhibited further enhanced expression levels in HMCLs in line with the trend in primary tumors, thus strengthening the putative involvement of these genes in sustaining the intensity/aggressiveness of the neoplastic phenotype.
A further finding of the present study was the identification of distinct transcriptional signatures in groups of pPCL that clearly showed different outcome. The 27 genes carrying a highly significant correlation with OS were prevalently upregulated (17/27, 63%) in the pPCL patients with the poorest outcome. The occurrence of this peculiar pattern is predictive of poorer OS, independent of the major cytogenetic alterations, hematologic parameters, and known multiple myeloma gene-risk models. Conversely, the low-risk signature is strongly associated with the ASCT procedure, which points to a more favorable prognosis overall. Further investigations on larger cohorts and validation sets are warranted to better elucidate whether a specific transcription pattern at diagnosis may be suggestive of the eligibility of the patients to ASCT with consequently associated better outcome. Interestingly, none of the 27 genes were recognized in known multiple myeloma high-risk signatures (24, 28, 30), whereas only a minor fraction were dispersedly included in the gene lists characterizing the 7 multiple myeloma molecular subgroups in the UAMS database (26). Interestingly, the platelet/endothelial cell adhesion molecule 1 (PECAM1), whose expression was most significantly associated with OS, encoded a member of the immunoglobulin superfamily involved in angiogenesis and activation of integrins. Notably, a recent work described its involvement in the lenalidomide effect in angiogenesis inhibition and the interaction with cadherin 5 and β-catenin, which is critical for endothelial cell cord formation (37). It could be speculated that the lower levels of PECAM1 associated with more aggressive disease might lead to less dependency of malignant PCs from microenvironment and reduced substrates for drugs, such as lenalidomide, that target microvessel formation.
The primary endpoint of the study was set to response after 4 cycle of therapy; 6 patients failed to reach response; the gene expression analysis highlighted 3 genes that clearly distinguished nonresponders from responders (albeit at different treatment response, from partial to complete) to the treatment with LD. Little is known about the biological role and the function of CYB5D2, EDEM3, and YIPF6, which leaves open the investigation of these genes as putatively involved in drug response or associated with PCs with a resistant phenotype. Additional studies on larger cohorts are warranted on this aspect.
Overall, our data suggest that pPCL might not represent a unique transcriptional and clinical entity, albeit clearly distinguishable from multiple myeloma on the basis of the expression of a large signature, and provide insights for further investigations of mechanisms underlying the biology of aggressive forms of PC dyscrasia.
Disclosure of Potential Conflicts of Interest
M. Offidani has received honoraria from the speakers bureau of Celgene and Janssen (Minor: $10,000 or less). No conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: L. Agnelli, A. Palumbo, A. Neri, P. Musto
Development of methodology: L. Agnelli, A. Palumbo
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K. Todoerti, S. Fabris, M. Lionetti, L. Mosca, V. Grieco, G. Bianchino, C. Mazzoccoli, L. De Luca, M.T. Petrucci, F. Di Raimondo, A. Falcone, P. Musto
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): K. Todoerti, L. Agnelli, M. Lionetti, G. Tuana, T. Statuto, F. Morabito, P. Omedè
Writing, review, and/or revision of the manuscript: K. Todoerti, L. Agnelli, F. Morabito, M. Offidani, F. Di Raimondo, P. Tassone, M. Boccadoro, P. Musto
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): V. Grieco, F. D'Auria, A. Neri
Study supervision: L. Lombardi, A. Palumbo, A. Neri
Grant Support
This work was financially supported by grants from: the Associazione Italiana Ricerca sul Cancro (AIRC) grant 5xmille #9980 2010-15 (to P. Tassone, A. Neri, and F. Morabito); the AIRC grant IG10136 and the Ministero Italiano dell'Istruzione, Università e Ricerca (MIUR) grant 2009PKMYA2 (to A. Neri); and the Fondazione Matarelli, Milan, Italy. M. Lionetti was supported by a fellowship from the Fondazione Italiana Ricerca sul Cancro (FIRC). P. Musto is supported by a Research fund from Celgene. K. Todoerti was supported in part by the Italian Health Minister, Finalized Research for Young Researchers, CUP Project E66110000230001.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.