High WEE1 expression is independently linked to poor survival in multiple myeloma

Introduction

Multiple myeloma (MM) is a hematologic malignancy associated with a malignant proliferation of plasma cells [1]. Although the disease is usually responsive to upfront therapies, MM remains incurable even in patients who achieve undetectable levels of disease, with relapse considered largely inevitable [2]. The genomic makeup of MM is highly heterogeneous, and different studies have identified multiple subtypes associated with varying prognostic outcomes using different data modalities [1, 3,4,5]. Standard methods to prognosticate the length of progression-free survival (PFS) include the International Staging System (ISS) [6], Revised ISS (R-ISS) [7], and the Second Revision of the ISS (R2-ISS) [8]. These tools rely on surrogates for disease burden and identification of specific tumor cytogenetic abnormalities. These scoring systems each have a PFS concordance index (c-index) below 60%, leaving room for improvement [9, 10].

In addition to providing genomic information, scoring systems informed by gene expression have been proposed for prognostication, including GEP70 and SKY92 [11, 12]. These expression-based signatures have shown potentially complementary information to ISS staging [11]. In [5], we conducted a large unsupervised genomic network study where we applied a novel measure of network connectivity, Ollivier-Ricci curvature, to RNA-sequencing (RNA-seq) and copy number alteration (CNA) data from newly diagnosed MM (NDMM) patients. We examined patterns of gene-gene interactions in MM and identified novel pathways and genes associated with poor prognosis. By examining the impact of gene expression via a network, we identified a novel eight-gene signature: BUB1, MCM6, NOSTRIN, PAM, RNF115, SNCAIP, SPRR2A, and WEE1. Of these eight genes, WEE1 was the only gene that was included in a previously published gene signature, GEP70 [13]. Furthermore, WEE1 was the most prognostic for PFS, suggesting it might play a role in MM. However, the role of the WEE1 in MM has not been thoroughly studied, and much remains unknown about its prognostic significance with respect to known biomarkers of MM.

WEE1 is a tyrosine kinase involved in multiple aspects of the cell cycle process, including the G1-S checkpoint, S phase, and G2-M checkpoint [14, 15], but believed to exert its most significant clinical impact in the G2-M checkpoint. For non-cancerous cells, DNA damage is often repaired at the G1-S checkpoint. In cancerous cells, the G1-S checkpoint may be deficient, and therefore, cancerous cells rely on the G2-M checkpoint for DNA damage repair [16]. In the G2-M checkpoint, WEE1 regulates cyclin-dependent kinase 1 (CDK1) [17,18,19], with high WEE1 expression suppressing CDK1 expression and maintaining the cell in a DNA repair state [20, 21]. Conversely, low WEE1 expression correlates with a rise in CDK1 expression, which allows the cell to enter mitosis [17].

WEE1 inhibition has been shown to dysregulate the cellular machinery associated with the first stage of mitosis in the G1-S transition [22] and can induce apoptosis by forcing mitotic entry [23]. For a cell to successfully complete the cell cycle, WEE1 expression levels must rise and fall in relation to each stage of the cycle. High WEE1 expression has recently been shown to be associated with disease aggressiveness in some solid tumors including breast cancer [24], ovarian cancer [25], and melanoma [26,27,28,29]. Several WEE1 inhibitors are currently in phase 2 clinical trials; these trials are evaluating the therapeutic efficacy of WEE1 inhibition [30, 31]. WEE1 inhibitors have also shown promise in other cancer types including sarcomas [32] and breast cancers [33], as well as hematological malignancies [34].

In MM, preclinical studies have shown promising results when inhibiting WEE1 in cell lines and mouse models in conjunction with other factors [35,36,37,38,39]. WEE1 inhibitors, in combination with bortezomib, can induce apoptosis in MM cell lines more efficiently than bortezomib alone [35, 36]. further, in [37], the authors show that bortezomib in combination with a DNA damage response (DDR) inhibitor targeting ATM/ATR/WEE1 triggers apoptosis. In [38], the authors examine the relationship between WEE1 and CHK1 in MM and report that targeting both kinases induces apoptosis in MM cell lines. In [39], the authors suggest targeting CTPS1 in conjunction with either CHEK1, ATR, or WEE1 inhibition can induce apoptosis in MM cell lines.

In this study, we show that high WEE1 expression defines a high-risk subtype of MM, independent of both known markers of MM and treatment types. WEE1 expression has comparable prognostic value as compared to the traditional MMISS. Additionally, high WEE1 expression is not reflected by corresponding changes in expression throughout the transcriptome. The high WEE1 expression subtype is characterized by dysregulation of the P53 pathway. Together, this work suggests that in a subpopulation of MM patients, WEE1 may play an outsized role and should be studied as a potential therapeutic target.

Methods

In this study, we applied a variety of bioinformatic and machine learning-based methods to MM datasets to examine the role of WEE1 in MM.

CoMMpass data

The RNA-seq and CNA data used is from the Multiple Myeloma Research Foundation’s CoMMpass dataset, release version 19. Further information on the data collection and curation methods has previously been published [40, 41]. The details of the patients selected for this study along with the preprocessing and feature computations are described in detail in [5]. Briefly, for inclusion in this study, subjects must have RNA-Seq and CNA data extracted from the bone marrow plasma cells before the start of treatment and both demographic and survival information available (N = 659). Gene inclusion was based on overlap with the Human Protein Reference Database [42].

Gene expression profiling (GEP) data

The GEP data used is from the University of Arkansas’s Total Therapy 2 (TT2, N = 341) and Total Therapy 3 (TT3, N = 214) trials. The details of these trials are described in [43, 44]. Briefly, the plasma cells were collected via a bone marrow biopsy of NDMM patients before treatment and gene expression profiling data were collected. TT2 and TT3 were different treatment regimens. Note that for this dataset, event-free survival was reported.

High-risk group membership

For each data modality—RNA-seq and GEP—patients’ WEE1 expression values were sorted in descending order and the top tertile was labeled as WEE1-high and the bottom tertile was labeled as WEE1-low. The center third was not considered in this study. The same process was independently repeated for CNA values.

Prognosis and confounder analysis

The prognosis was modeled using Kaplan–Meier (KM) survival curves for PFS. To determine the effect of WEE1 relative to known biomarkers of MM, we used a multivariate Cox proportional hazards model [45] with the RNA-Seq data to predict PFS. In it, we modeled nine markers: hyper APOBEC, chromothripsis, hyperdiploidy, MAF translocation, MYC translocation, t(4;14), t(11;14), TP53 mutation, and gain 1q21. As outlined in [5], hyperdiploidy was defined by more than 2 gains involving >60% of the chromosome affecting chromosomes 3, 5, 7, 9, 11, 15, 19, or 21. Mutational signatures were assessed using mmsig (https://github.com/UM-Myeloma-Genomics/mmsig), a fitting algorithm designed for MM to estimate the contribution of each mutational signature in each sample [46]. APOBEC-mutational activity was calculated by combining SBS2 and SBS13, with the top 10% being defined as hyper APOBEC [47, 48]. The complex structural variant chromothripsis was defined by manual curation according to previously published criteria [49]. High-risk and low-risk groups were analyzed separately to see which factors differed between the groups. To show the prognostic effect of WEE1, irrespective of known biomarkers, KM survival curves for PFS stratified by each factor were plotted.

Machine learning analysis

We used random survival forests [50] to determine the prognostic value of WEE1, its gene network neighbors, and ISS. Briefly, random survival forests offer the advantages of random forests with the addition of incorporating survival information including event duration and censorship information. WEE1 neighbors were extracted from the STRING database [51]. WEE1 neighbors were defined as genes which have a known interaction with WEE1 with a probability greater than 0.7. The neighboring genes were considered to see if changes in WEE1 expression were reflected by changes in the expression of known interacting genes. ISS staging was provided by the CoMMpass dataset. We used the concordance index (c-index) as the evaluation metric. WEE1 expression was predicted using random forest regression models to see if neighboring genes contained signals relevant to the abnormal increase in WEE1 expression. Feature importances were computed using the permutation importance method in sci-kit-learn and the 15 most importances are reported [52]. The full parameter details of the models used are available on GitHub (www.github.com/aksimhal/WEE1-in-MM). Models were evaluated using five-fold cross-validation repeated ten times.

Differential gene expression analysis

To see differences in patterns of gene expression between the WEE1-high and WEE1-low cohorts, we computed the differential gene expression using DESeq2 [53]. The p-values from this analysis were corrected for multiple hypothesis testing using the BH-FDR method. Genes with a corrected p value less than 0.05 and an absolute log2 fold change greater than two were considered significant. To see which pathways become dysregulated in WEE1-high, we used the Gene Set Enrichment Analysis tool to evaluate the selected genes [54, 55]. The utilized pathways are from the hallmark gene set collection from the human molecular signatures database [56].

Ethics statement

This project was performed in accordance with the Declaration of Helsinki and the Belmont Report. All data referenced have been previously published and made publicly accessible. MMRF CoMMpass data was anonymized prior to access. Patient enrollment and data safety were managed by the MMRF per their organization guidelines. At no point did any investigator have access to protected health information for any included patients. Anonymized GEP datasets were accessed through the National Institutes of Health’s Gene Expression Omnibus public-facing website.

Results

Data overview

Genomic and clinical characterization of MM outcomes were stratified by WEE1 expression using the CoMMpass dataset (N = 659). The mean age was 62.5 ± 10.7 years and 60% were male; ISS distribution was 35/35/30%, and 53% received an autologous stem cell transplant (ASCT). An overview of the differences between the WEE1-high and WEE1-low groups is provided in Table 1. While some of the known markers of MM are significant between the two groups, including age, hyperdiploidy, t(11;14), MAF and MYC translocations, chromothripsis, hyper APOBEC, gain 1q21, and TP53 mutational status, ISS is not. In the WEE1-high cohort, 25% were R-ISS I, 62% were R-ISS II, and 13% were R-ISS III. In the WEE1-low cohort, 29% were R-ISS I, 64% were R-ISS II, and 7% were R-ISS III. For the validation datasets, TT2 and TT3, baseline clinical data and gene expression data were available. For TT2, the mean age was 56.3 ± 9.8 years and 57% male; for TT3, the mean age was 58.6 ± 8.8 years and 67% male.

Table 1 Difference in CoMMpass data patient characteristics between WEE1-high and WEE1-low cohorts.

Full size table

WEE1 is prognostic for outcomes in RNA-seq and GEP datasets

In the RNA-seq data from the CoMMpass dataset, differences in PFS between WEE1-high and WEE1-low cohorts are statistically significant (p < 1e-9), as shown in Fig. 1A. These results are validated in the TT2 and TT3 datasets (Fig. 1B, C). Note, when stratifying subjects by the top tertile and bottom tertile of CNA data from the CoMMpass dataset, the effect is not observed.

**Fig. 1: Prognostic value of *WEE1* expression from RNA-seq and GEP data.**

Multivariate modeling shows that WEE1 is an independent prognostic factor in MM

Multivariate Cox proportional hazards modeling shows that the prognostic effect of WEE1 is independent of known MM markers, including those shown to be significant in Table 1. The prognostic effect is independent of hyperdiploidy, t(4;14), t(11,14), TP53 status, as well as emerging risk factors, the complex structural variant chromothripsis and APOBEC-mutational activity, shown in Fig. 2A and Supplementary Table 1A. When examining only the WEE1-high cohort, none of the markers significantly predicted PFS (Fig. 2C, Supplementary Table 1B). Similarly, in the WEE1-low, none of the markers significantly predicted PFS (Fig. 2B, Supplementary Table 1C).

**Fig. 2: Cox proportional hazards (CPH) modeling of MM markers and *WEE1* expression.**

WEE1 is prognostic for outcomes independent of known biomarkers

The WEE1-high and WEE1-low cohorts have statistically significantly different PFS outcomes when stratifying for each known MM marker. KM plots show significant separation when looking at groups defined by the presence of hyperdiploidy, t(11;14), MAF and MYC translocations, chromothripsis, and TP53 deletion (Fig. 3, Supplementary Fig. 1). KM plots were also significant when looking at the groups defined by the lack of a known MM marker (Fig. 3, Supplementary Fig. 2). WEE1 cohort membership differentiates PFS outcomes by an average of 1.98 years in cohorts with a marker, and 2.18 years in cohorts without the marker (Table 2).

**Fig. 3: Kaplan–Meier curves stratified by MM markers show the prognostic signal in *WEE1* expression.**

Table 2 The difference in median progression-free survival (PFS) is based on a given biomarker.

Full size table

WEE1 is prognostic for outcomes independent of treatment type

The WEE1-high and WEE1-low cohorts have statistically significantly different PFS outcomes when stratifying the treatment options listed in the CoMMpass dataset. ASCT, bortezomib/immunomodulatory agents (IMIDs), bortezomib, and carfilzomib/IMIDs cohorts were all significantly different when stratified by WEE1-high and WEE1-low (Fig. 4). The mean difference in PFS is 1.91 years.

**Fig. 4: Kaplan–Meier curves show the effect of *WEE1* expression on treatment type.**

WEE1 expression has comparable prognostic value as ISS

RNA-seq-based WEE1 expression has a comparable prognostic value (c-index: 0.58 ± 0.04) as ISS (c-index: 0.61 ± 0.03). Combining WEE1 and ISS has a c-index of 0.63 ± 0.03.

WEE1-high cohort is 3.2× less predictable than the WEE1-low cohort

As WEE1 expression increases, the relationship between WEE1 and genes known to interact with WEE1 becomes dysregulated. When modeling WEE1 expression with known interacting genes, the prediction error increases by 3.2 times between the WEE1-high and WEE1-low cohorts. In the WEE1-low cohort, the known interacting genes that contribute more than 5% to the prediction are CDK1, CHEK1, CDT1, AURKB, and PLK1 (Fig. 5A). In the WEE1-high cohort, the genes are CDC25B, HSP90AA1, CDK6, PLK1, CDR2, SKP2, and CDK2 (Fig. 5B).

**Fig. 5: Random forest feature importance plots.**

P53 pathway-related genes are differentially expressed between WEE1-high and WEE1-low cohorts

A differential gene expression analysis between the WEE1-high and WEE1-low groups identified 146 overexpressed genes and five underexpressed genes. Overexpressed genes are part of three pathways: P53, downregulated UV response, and mitotic spindle. Only five genes were underexpressed: FPR1, IFNA5, LRP2, POU2F3, and RAB11FIP1. The full list of differentially expressed genes is in Supplementary Table 2 and illustrated in Supplementary Fig. 3.

Discussion

Prognostic markers in MM rely on either assessment of tumor burden or specific cytogenetic abnormalities; transcriptional characteristics of myeloma are not currently considered in this setting. Here, we have identified that high WEE1 expression represents an independent biomarker prognostic of poor outcomes in NDMM, and that this effect is independent of known cytogenetic risk factors and treatment strategies (Figs. 1–3). This includes the common metric of staging—ISS. Even in cohorts defined by a MM risk factor, WEE1 expression contains an additional prognostic signal that further differentiates outcomes. Random survival forest modeling showed that WEE1 expression alone has as much prognostic power as ISS staging. These findings were seen both retrospectively using the CoMMpass dataset and independently validated in two additional MM datasets. Differential gene expression analysis showed that the P53 pathway is the most significantly affected pathway in the WEE1-high cohort. Furthermore, the outcomes stratification effect described here is only observed when stratifying by gene expression, not by gene CNA.

Random forest modeling of the local WEE1 genomic network showed that the overexpression of WEE1 is not correlated with an increase or decrease in any genes locally connected with WEE1. Increased WEE1 expression was not reflected in a rise in the expression of any other cell cycle kinases, such as PLK1 or CDK1. Random forest modeling of the low-risk group showed an association with CDK1, which follows known biology. However, in our defined high-risk group, CDK1 was not in the top 15 genes most associated with the high-risk WEE1 signal. This further suggests that WEE1 expression represents an independent prognostic marker that is likely not merely reporting on another known cytogenetic risk factor.

WEE1 is a key player during the cell cycle, and its specific roles in the S phase and the G2-M checkpoint are well documented. WEE1 acts as a tumor suppressor gene in certain types of breast cancer. However, for the majority of solid and blood cancers, such as ovarian cancer and acute lymphoblastic leukemia, WEE1 acts as an oncogene. Further work is needed to understand the role of increases in WEE1 expression in MM as these findings can enable new WEE1-directed treatments in MM patients with MM and other malignancies.

Of note, differences in PFS among patients with TP53 deletions when stratifying by WEE1 expression were remarkably large. Patients with TP53 deletions often have the poorest clinical outcomes with MM treatment across multiple published datasets. Additionally, differential gene expression analysis between the high-risk and low-risk groups showed that genes associated with the hallmark P53 pathway were differentially expressed. TP53 regulates DNA damage in the G1-S checkpoint. Faulty P53 function may lead to a larger reliance on WEE1 activity to maintain genomic integrity. If both TP53 and WEE1 are abnormal, it is possible that DNA repair becomes dysfunctional.

We have demonstrated that stratification of MM patients with TP53 deletions by MM cell WEE1 expression may represent an alternative method of risk-stratifying patients. Additionally, our data suggests that WEE1 inhibition may be especially effective in patients with altered P53 pathways, though further investigation is needed to identify if the observed association is causal. There are currently five WEE1 inhibitors in clinical trials [57] for other cancer types which will advance our understanding of the efficacy of WEE1 inhibition, the exact mechanism of its actions, as well as a possible new treatment option for MM patients.

Our results show that WEE1 expression is prognostic independent of known biomarkers, differentiates outcomes associated with known markers, is upregulated independently of its interacting neighbors, and is associated with dysregulated P53 pathways. This suggests that WEE1 expression levels may have clinical utility in prognosticating outcomes in NDMM and may support the application of WEE1 inhibitors to MM preclinical models. Determining the causes of abnormal WEE1 expression may uncover novel therapeutic pathways.