Progression signature underlies clonal evolution and dissemination of multiple myeloma

Cancer is a genetically complex evolutionary process whereby transformed cells continue to acquire genetic and/or epigenetic lesions, giving rise to heterogeneous populations of functionally distinct clones.^1-3  The number of mutant clones per cancer can range from tens to hundreds (most commonly) or thousands, out of which only few clones are functionally relevant true drivers.^1,3  Forced to compete for limited resources under natural selection, with the presence of microenvironmental and other constraints, the most resilient clones with best fitness to self-renew are selected for extensive proliferation, migration, and invasion.⁴  Ultimately, cancer cells that survive can go on to repopulate distant tumor microenvironments, possibly carrying new alterations that enhance their malignant potential.^3,5,6 

Multiple myeloma (MM), an incurable plasma cell malignancy of the bone marrow (BM), is a particularly heterogeneous type of cancer^7-10  whose progression from the well-defined precursor stages of monoclonal gammopathy of undetermined significance (MGUS) and smoldering MM (SMM) is underlain by marked evolution.^11-13  As such, although at the early stage of disease, a few clones with limited number of alterations are present, subsequent abnormalities are acquired as patients progress to overt disease, conferring fitness advantage and allowing for tumor expansion.^8,14  Nevertheless, genetic/epigenetic aberrations are seldom enough to drive cancer progression by themselves; a permissive microenvironment and reciprocal interactions between tumor and microenvironmental cells are equally important.^15,16  Multiple components in the tissue environment can influence cancer clonal evolution, and cancer cells in turn can also remodel the microenvironment for their competitive advantage.^15,17,18  Thus, a model of MM progression with the presence of tumor microenvironment is crucial and can provide a unique platform to understand clonal heterogeneity, evolution, and dissemination mechanisms.

In this study, we describe a novel “bone chip” xenograft mouse model (MM PrEDiCT), which in combination with a fluorescent or DNA-barcodes system can be used to track and profile heterogeneous clonal evolution as they disseminate from the primary sites to peripheral blood and distant BM sites. Subsequently, we use this model to define an RNA expression signature of progression, which we find to correlate with disease progression and survival in MM patient data, confirming its relevance. Finally, we identify potential drivers of progression, HMGA1 and PA2G4, in our model and validate their significance in MM disease progression and dissemination through an in vivo CRISPR screen and in vitro assays.

All mouse experiments were performed with an Institutional Animal Care and Use Committee–approved animal protocol at our facility (Dana-Farber Cancer Institute). Severe combined immunodeficiency (SCID)-beige mice (C.B-Igh-1b/GbmsTac-Prkdc^scid-Lyst^bg N7) were obtained from Taconic. For tumor-bearing bone chip implantation, femurs were resected from 6- to 7-week-old female SCID-beige mice. Two million myeloma cells were injected into the BM cavity of femurs resected from donor mice and were then transplanted subcutaneously into syngeneic recipient mice of the same age.

For RNA sequencing, poly(A) selection and complementary DNA synthesis were performed, followed by library preparation using Illumina TruSeq RNA Sample Prep Kit, sequencing (75-bp paired reads), and sample identification for quality control. RNA-sequencing data were processed by Kallisto,¹⁹  a pseudoalignment tool used to quantify abundances of transcripts, which were represented as transcripts per million. Estimated counts were imported by R package tximport for downstream analysis. Differential expression analysis was performed by DESeq2²⁰  with a false discovery rate (FDR) cutoff of 5%.

All mice were treated, monitored, and sacrificed in accordance with an approved protocol of the Dana-Farber Cancer Institute Animal Care and Use Committee.

BM samples from relapsed/refractory MM patients were obtained under the approval from the Dana-Farber Cancer Institute Institutional Review Board. Informed consent was obtained from all patients in accordance with the Declaration of Helsinki.

All statistical analysis was performed in R (https://www.r-project.org/), and multiple hypothesis testing was corrected using the Benjamini-Hochberg method.

To investigate the dynamics of clonal heterogeneity in vivo within the BM microenvironment, we developed a tumor dissemination xenograft model characterized by subcutaneous implantation of myeloma-bearing femur bone chips into syngeneic SCID-beige mice (supplemental Figure 1A, available on the Blood Web site). Two weeks following bone chip implantation, engraftment was assessed using confocal microscopy. At that point, functional vessels and bone remodeling could be observed, while myeloma cells had started colonizing the proximal vessel plexus, suggesting that the implanted bone chips had rapidly engrafted and were functioning as a permissive microenvironment for the growth and migration of myeloma cells (supplemental Figure 1B). Disseminated MM cells had colonized distant BM sites by week 6 (supplemental Figure 1C), while by weeks 8 to 10, host mice presented with limb paralysis, displaying multiple skeletal lytic lesions. We were able to reproducibly model distant BM dissemination of various myeloma cell lines using bone chip implants (supplemental Figure 1D), mimicking the human MM dissemination phenotype. In contrast, direct subcutaneous injection of MM.1S into SCID-beige mice did not reveal significant amount of disseminated GFP⁺ cells in the blood, BM, or spleen (supplemental Figure 2). These data suggest that MM cells alone are not sufficient for progression, and a permissive BM stromal environment is required for disease progression/dissemination. Thus, we named this model MM PrEDiCT (MM progression through evolution and dissemination of clonal tumor cells).

Subsequently, we took advantage of the fluorescent tagging system to generate a stable 15-color cell library with cells bearing fluorescence proteins encoding GFP, BFP, RFP, and iRFP (Figure 1A).²¹  These 15-color subpopulations can be flow-sorted, mixed in equal proportions (Figure 1), and “visually” tracked in the PrEDiCT model, hereafter referred to as the “rainbow” system.

Rainbow cell–bearing bone chips were implanted as described above. Approximately 60 days after transplantation, rainbow cells from in vitro culture, in vivo primary bone chips, and distant BM were isolated and analyzed by flow cytometry. Results revealed that both primary and distant tissues are composed of persistent and fluctuating clones. The fluorescence distribution of in vitro culture was approximately even at days 30 and 60 after transplantation (Figure 1B). By contrast, although all 15 subpopulations in the primary tumor could be detected at day 60, they showed uneven distributions in size, with some colors unequivocally overtaking others (Figure 1C), reflecting the existence of clonal competition even within the primary bone chip, where few “dominant” clones adapted quickly and were able to evolve and outcompete the rest. Upon disease progression, reduced clonal (color) diversity was observed in the distant BM and a few large clones, accounting for majority of all cells (Figure 1C; supplemental Figure 3A). Interestingly, clones that dominate the matched circulating tumor cell (CTC) population were also found to dominate distant BM sites (Figure 1C), consistent with certain clones locally invading the primary tumor parenchyma, intravasating into nearby microvessels, and seeding distant BM sites. Strikingly, although “dominant” clones exhibited intermouse variability, similar “dominant” clones were observed in the left and right femur of each mouse, with tightly correlated fluorescence distributions between matched femurs and CTCs (supplemental Figure 3B). These results suggest that following circulatory dissemination, the tumor microenvironment might be largely similar across different BM sites, allowing for expansion of similar clones with minimal site heterogeneity.

We next sought to understand the clonal architecture within each site by observing the implanted bone chip and distant BM sites with confocal microscopy. Tumor cells appeared to form multiple clusters in both sites. Interestingly, under high magnification, the primary site clusters appeared to comprise various cell subpopulations, as evidenced by the multitude of colors seen in those clusters (Figure 1D), whereas most of distant BM clusters comprised single-color populations (Figure 1E). These results demonstrate the existence of significant spatial heterogeneity within each distant BM site, with clear demarcation of each clone’s area of residence.

To improve clone-tracking depth and accuracy, we developed a lentiviral DNA-barcoding system with a 12-base random sequence tag that can generate ≤16 million unique and inheritable barcodes detectable by next-generation sequencing.²²  Two million barcoded cells were injected into femurs as described above. As expected, while the pretransplantation pool showed no major clonal bias (Figure 2B), primary tumors, disseminated cells, and CTCs exhibited clonal selection, with distant sites demonstrating much less diversity (Figure 2C). In addition, within each mouse, left/right femurs and CTCs exhibited the same “dominant” clones (Figure 2C), while clonal distributions in different mice were distinct, as evidenced by strong shifts in the respective cumulative distribution rankings (Figure 2D). These results further confirm that clonal competition, in the context of tumor heterogeneity, started at an early stage, and only a small number of clones that can adapt to the environment and evolve quickly to compete with other clones have the ability to successfully disseminate.

We also compared different injection methods on clonal evolution. The clonal distributions of CTCs and disseminated BMs of both intravenous and intratibial injections had much higher diversity (supplemental Figure 4), indicating that clonal selection/competition power was much lower than bone chip model. The implanted bone chip acts as a strong bottleneck that selects for cells that are able to first survive at the primary sites and it is out of those that circulation is permitted, enabling a clearer tracking of winning clones from primary tumor to circulation and to distant BM. Thus, PrEDiCT model largely reflects a multistage process of dissemination/metastasis, including invading the local primary tissue, intravasating into the tumor vasculature, circulating, and ultimately extravasating into the parenchymal of distant tissue for further colonization.^23-25 

To explore the molecular mechanisms underlying disease progression, we performed RNA sequencing on 3 human cell lines (MM.1S, IM-9, and OPM2) harvested from matched primary bone chips and distant BMs separately. Since IM-9 is an Epstein-Barr virus–transformed B lymphoblastoid cell line derived from MM patients and MM.1S and OPM2 have notable genetic differences, such as mutation status of the TP53 gene, we thus used MM.1S as our main model for signature discovery and OPM2 and IM-9 as references. Differential expression analysis identified 1109 upregulated genes and 1865 downregulated genes in the MM.1S model at an FDR cutoff of 5% (supplemental Table 1). The top 300 up- and downregulated genes were selected as a signature to track the dissemination potential of MM cells (Figure 3A), and we designated this gene set as the progression signature. We employed gene set enrichment analysis and showed that both up- and downregulated signatures were significantly enriched in the IM-9 model (supplemental Figure 5A-B), and only downregulated signatures were also significantly repressed in the OPM2 model (supplemental Figure 5C).

Next, we sought to validate the relevance of this signature in human MM progression using a public patient gene expression data set (GSE6477).²⁶  Strikingly, our signature was sequentially activated during MM progression, with more significant P values observed at later stages (Figure 3B). Of note, signature enrichment could already be detected at the MGUS stage (P = .004), potentially reflecting the malignant potential of this early precursor stage, and was further increased in plasma cell leukemia (Figure 3C), an advanced and aggressive form of MM. Most importantly, the progression signature was significantly associated with inferior overall survival (P = .0344) (GSE24080; Figure 3D; supplemental Figure 6).²⁷  Together, our results suggest that the PrEDiCT model closely mirrors MM disease progression within the patient clinical setting at the phenotypic, cellular, and molecular level. Pathway enrichment analysis of dysregulated genes in MM.1S indicates that several pathways contribute to progression potential of MM cells (Figure 3E; supplemental Table 2). For example, apoptosis and the p53 signaling pathway were repressed, which are related to proliferation and transformation. Some immune pathways, such as Fc γ receptor−mediated phagocytosis signaling, are repressed as well. Spliceosome is the only activated pathway, and its biological significance remains to be explored.

Gene expression signatures are known to represent phenotypical markers and do not usually suggest upstream regulators. Previous studies have shown that regulon analysis using the master regulators (MRs) inference algorithm can help identify aberrantly activated tumor drivers.²⁸  Briefly, we applied ARACNe²⁹  on a public MM patient expression data set (GSE6477) and identified a regulatory network, which contains 591 MRs targeting 12 395 genes, with a P value cutoff of 10⁻⁸. VIPER²⁸  was then applied to identify potential MRs that could be driving MM progression and dissemination in our model, as measured by enrichment of progression signature. A total of 28 MRs were found to be significant with an FDR threshold of 5%, including 15 activated and 13 repressed regulators (Figure 4A; supplemental Table 3). Of note, some of the repressed MRs were downregulated in distant BM sites, including FOSB, JUN, and KLF6. Next, we looked at the MR gene dependency scores in the Project Achilles data set (v20Q2) and found that the majority of the MRs discovered in our model are indeed essential in MM cell lines (Figure 4B; supplemental Figure 7).³⁰  For the purposes of validation, we selected activated ones, as they may serve as novel therapeutic targets. After ruling out previously reported factors in MM like MYC,^31-33 CDKN2A,^34-36 TOP2A,^37,38  and TCF3,³⁹  as well as nonessential genes in MM based on DepMap data, HMGA1, PA2G4, and TRIM28 were selected based on statistical significance, novelty, and putative functional importance as potential therapeutic targets (Figure 4C).

We designed a targeted CRISPR library containing 30 single guide RNAs (sgRNAs) targeting HMGA1, PA2G4, and TRIM28, as well as 100 nontargeting control sgRNAs (supplemental Table 4). Lentiviruses expressing the library were used to transduce MM.1S-GFP-Luc-Cas9 cells, which were then cultured in vitro for 1 week before being transplanted into SCID-beige mice. For each gene, we sequenced and compared distribution of normalized sgRNA counts in both primary tumors and BM at the time of limb paralysis (10-12 weeks after transplantation) using the software package MAGeCK (Model-Based Analysis of Genome-wide CRISPR/Cas9 Knockout).⁴⁰  sgRNAs targeting HMGA1 and PA2G4 were first significantly depleted in the late-stage primary sites compared with injection pool (supplemental Figure 8A) and further depleted in distant BM sites, as compared with primary sites, while nonspecific sgRNAs have similar median read counts in primary tumor and BM (Figure 4D). Both HMGA1 and PA2G4 were upregulated MRs, which tend to be essential for cell survival in MM cell lines (Figure 4B), indicating both genes can be potential targets for treating myeloma. When looking at the gene expression data from different stages of MM using GSE6477 data set, expression levels of both HMGA1 and PA2G4 were found to be significantly higher in newly diagnosed and relapsed MM (P = .011 and P = .013, P = 3.90e-05, and P = 1.23e-04, respectively) compared with healthy donors, and in the case of PA2G4, expression levels were even significantly higher in MGUS and SMM (P = .019 and P = .04, respectively) (Figure 5). Our results suggest that these genes become increasingly important to MM cells as the disease progresses, confirming their importance in myeloma progression and as potential therapeutic targets.

We next tested whether HMGA1 and PA2G4 regulate the same set of genes from progression signature. Loss of HMGA1 and PA2G4 were confirmed by western blot (supplemental Figure 8B). Real-time polymerase chain reaction was performed on these cells for 5 significantly upregulated (CD38, FANCF, MGST1, AARSD1, and LBR) and 5 downregulated (KLF6, LGMN, SIRT2, GAB1, and TIPARP) genes shared by the MM.1S, OPM2, and IM-9 models. Interestingly, among those genes activated at disseminated sites, CD38, FANCF, and LBR had significantly reduced expression in HMGA1-deficient cells, but not in PA2G4-deficient cells (supplemental Figure 8C). Expression levels of both MGST1 and AARSD1 were slightly increased in both HMGA1- and PA2G4-knockout cells. Similarly, for the 5 genes downregulated at disseminated BM sites, loss of HMGA1 and PA2G4 drastically increased the expression of LGMN, SIRT2, and TIPARP, but only loss of PA2G4 increased KLF6 expression (supplemental Figure 8D). These data suggest that these 2 MRs potentially regulate different gene networks to give cells distinct advantages for disease progression; in other words tumor heterogeneity gives rise to subclones that possess different evolutionary advantages over the others within the tumor environment.

To further confirm the important roles of HMGA1 and PA2G4 in driving MM progression, we performed in vitro proliferation, migration and adhesion assays with and without BM stromal cells (BMSCs) isolated from relapsed/refractory MM patients. The overall proliferation rates of MM.1S were increased when cocultured with BMSCs; however, loss of HMGA1, PA2G4, and MYC caused significant reductions in proliferation when compared with control cells in both culturing conditions (Figure 6A). We used MYC here as a positive control. Similarly, all knockout cells all showed significantly reduced migration and adhesion (Figure 6B and 6C) toward BMSCs, confirming the important roles of these genes in driving MM progression and dissemination. Interleukin-6 is a proinflammatory cytokine crucial to the growth, proliferation, and survival of myeloma cells.^41-43  Addition of 10 ng/mL interleukin-6 increased overall proliferation and migration; however, loss of HMGA1, PA2G4, and MYC compromised the proliferation and migration ability of MM cells significantly (supplemental Figure 9).

MM is an incurable plasma cell malignancy of the BM. Although many alterations regulating MM disease progression are tumor cell autonomous, be they genetic or epigenetic,^8,44-46  they are not the sole determinants of tumor behavior and are insufficient to induce progression and dissemination.^47,48  A permissive microenvironment is required for overt malignancy to emerge.^47,49  Indeed, studies have shown that the tumor microenvironment is a key regulator in many steps of the invasion-metastasis cascade, including oncogenesis, egress, protection in the circulation, preparation of the metastatic niche, organ-specific homing, and tumor colonization.⁵⁰  As such, disease progression in MM is the result of a well-orchestrated, complex network of interactions that cannot be studied outside of the context of the BM microenvironment. Our goal in this study was to develop a xenograft mouse model to reproduce the clonal evolution underlying dissemination/metastasis in progressing patients in the context of stromal microenvironment and enable the investigation of mechanisms underlying these events where mathematical analysis and library or drug screening can be further performed.

Studies by many groups have made the effort to track clonal outgrowth and dissemination in syngeneic mouse models using intravenous injection method, and similar conclusions were made that establishment of dissemination/metastasis is extremely inefficient and individual clones exhibited marked differences in proliferation and gene expression patterns in BM niche.^51-53  We also compared the differences in clonal evolution using bone chip, intravenous, and intratibial injection methods and observed significantly reduced clonal diversity in CTCs and distant BM in the PrEDiCT model, indicating that clonal selection/competition power was much higher. The implanted bone chip acts an optimal niche but a strong bottleneck that selects for cells that first survive at primary sites, and it is out of those that circulation is permitted. PrEDiCT could model the entire dissemination/metastatic cascade (invasion, intravasation, circulation, extravasation, and colonization). This approach revealed 3 major findings: (1) clonal competition was already present at the primary sites; (2) establishment of CTCs was extremely inefficient, and the efficiency was further decreased at distant BM sites; and (3) the individual disseminated MM cells exhibited marked differences in their proliferative fates, with the majority of final tumor burden within a bone being attributable to the progeny of very few clones. Thus, the PrEDiCT model provides a powerful platform for studying the role of tumor cell–autonomous alterations and BM stromal compartment in MM disease progression.

In an effort to elucidate expression changes associated with progression, we performed RNA sequencing on tumor cells (MM.1S, IM-9, and OPM2) isolated from primary and distant BM sites at the time of limb paralysis. We identified significantly up and downregulated genes shared across cell lines and designated those found significant in MM.1S cells as “progression signature.” In light of the data that tumor cells alone are not sufficient for dissemination/metastasis without the presence of a BM environment, the progression signature from the PrEDiCT model would reflect the changes in adaptation to the primary and secondary bone environment as well as the intrinsic evolutionary advantages evolved to enable dissemination during progression. Confirming our model’s relevance in studying human MM progression, gene set enrichment analysis of patient expression data showed that the progression signature was sequentially activated during MM progression, and more significant P values were observed when later stages of MM were compared with normal controls. Even more importantly, our signature was shown to be significantly associated with inferior overall survival in patients. These results suggest that the PrEDiCT model adequately represents the biological processes underlying disease progression in MM and thus can be used as a tool to improve our understanding of the mechanisms involved and identify new therapeutic targets.

To demonstrate the potential of our model for therapeutic target detection, we employed a MR analysis approach to identify potential upstream regulators of the progression signature. We selected and validated activated MRs, since they might serve as potential therapeutic targets. Loss of HMGA1 and PA2G4 compromised tumor adaptation to the primary sites and dissemination to distant sites in vivo, and in vitro assays further confirmed compromised abilities of these cells to proliferate, migrate, and adhere to BMSCs, indicating a role for HMGA1 and PA2G4, as progression-related MM additions and loss of both genes compromised the overall ability of tumor cells to overcome each bottleneck in metastasis process. Although HMGA1 has been implicated in oncogenic transformation, invasion and metastasis in various cancers,^54-60  its role in pathogenesis and progression of MM is largely unknown. Similarly, PA2G4 has been reported as a regulator for proliferation, but its role in cancers has being paradoxical.^61-64  Here, we report for the first time that HMGA1 and PA2G4 were sequentially activated along human MM progression in a public data set⁶⁵  and might play key roles in regulating dissemination of MM tumor cells within BM microenvironment, potentially as markers and therapeutic targets for progression.

Although both HMGA1 and PA2G4 were validated to be important for MM progression, the gene expression networks regulated by them are very different. Loss of HMGA1, but not loss of PA2G4, caused reduced expression of CD38, FANCF, and LBR, which are activated genes identified from the progression signature, and both genes could only inhibit some of the repressed genes (LGMN, SIRT2, and TIPARP) in the progression signature. This is in line with our clonal tracking observations that tumor clones at metastatic sites had much reduced diversity but still exhibited marked differences, with some clones dominating disseminated sites. In other words, we observed significant tumor heterogeneity from primary tumors to the metastatic sites, which gives rise to subclones that possess an evolutionary advantage over the others within the tumor environment, as manifested by distinct gene expression patterns, metabolism, motility, proliferation, and metastatic potential.

Recent studies showed that tumors are more than insular masses of proliferating cancer cells.^1,66,67  Instead, they are complex tissues composed of multiple distinct cell types participating in heterotypic interactions with one another.^4,46,68-70  One of the limitations of current model is lacking an intact immune system, while it is technically and experimentally challenging to perform bone chip modeling in syngeneic immune competent mice, we found a permissive BM stromal environment is crucial, because without it, dissemination is not achievable. Indeed, coculture of patient BMSCs with MM cells increased the malignant potential of MM cells, as manifested by significantly enhanced proliferation, migration, and adhesion. It will be interesting to further explore the specific components of stromal compartment and their biological functions in driving dissemination to accelerate our understandings of the complex BM environment in the future.

In summary, we established a novel xenograft mouse model termed PrEDiCT to track clonal dynamics and dissemination of MM, and we further validated HMGA1 and PA2G4 as potential drivers. Our model is unique in its ability to largely mirror the phenotypic, cellular, and mechanistic features of myeloma progression in the context of the BM stromal microenvironment. It constitutes a powerful platform for studying tumor cell and microenvironmental alterations in myeloma progression while providing new opportunities for investigation of mechanisms underlying these events where therapeutic target discovery or drug screening can be further performed.

This work is supported by National Institutes of Health, National Cancer Institute grants R01 CA181683-01A1 and R01 CA205954-01, Leukemia and Lymphoma Society grant 7004-14, and the Multiple Myeloma Research Foundation. A.M.R. was supported by Fondazione Associazione Italiana per la Ricerca sul Cancro (AIRC), Fondazione Regionale per la Ricerca Biomedica (Regione Lombardia), and JTC 2016 ERA-NET Transcan-2. Y.J.S. received a Young Investigator Award from the International Myeloma Society.

Contribution: Y.J.S., Y.M., J.S., and I.M.G. designed the study and wrote manuscript; Y.J.S., Y.M., J.S., M.M., S.M., A.M.R., A.S., Y.-T.T., Y.K., N.K.S., and B.B. performed the experiments; and Y.J.S., Y.M., J.S., R.S.-P., R.A.R, J.G.D., D.E.R., F. Mercier, F. Michor, D.T.S., and I.M.G. analyzed and interpreted the data.

Conflict-of-interest disclosure: I.M.G. has a consulting/advisory role with GlaxoSmithKline, Sanofi, Janssen, Celgene, Takeda, Karyopharm, AbbVie, GNS, Cellectar, Medscape, Genetech, Adaptive, Bristol-Myers Squibb, Aptitude, Curio Science, Magenta, and Oncopeptides. A.M.R. is on the advisory board for Amgen and Janssen Pharmaceuticals and has received research funding from AstraZeneca. J.G.D. consults for Tango Therapeutics. D.T.S. is shareholder for Magenta Therapeutics and serves as a director of Agios Pharmaceuticals, Magenta Therapeutics, and Editas Medicines. The remaining authors declare no competing financial interests.

Correspondence: Irene M. Ghobrial, Dana-Farber Cancer Institute, 450 Brookline Ave, Boston, MA USA; e-mail: irene_ghobrial@dfci.harvard.edu.