ERBB1/EGFR and JAK3 Tyrosine Kinases as Potential Therapeutic Targets in High-Risk Multiple Myeloma

Our main objective was to identify abundantly expressed tyrosine kinases in multiple myeloma (MM) as potential therapeutic targets. We first compared the transcriptomes of malignant plasma cells from newly diagnosed MM patients who were risk-categorized based on the patient-specific EMC-92/SKY-92 gene expression signature values vs. normal plasma cells from healthy volunteers using archived datasets from the HOVON65/GMMG-HD4 randomized Phase 3 study evaluating the clinical efficacy of bortezomib induction/maintenance versus classic cytotoxic drugs and thalidomide maintenance. In particular, ERBB1/EGFR was significantly overexpressed in MM cells in comparison to normal control plasma cells, and it was differentially overexpressed in MM cells from high-risk patients. Amplified expression of EGFR/ERBB1 mRNA in MM cells was positively correlated with increased expression levels of mRNAs for several DNA binding proteins and transcription factors with known upregulating activity on EGFR/ERBB1 gene expression. MM patients with the highest ERBB1/EGFR expression level had significantly shorter PFS and OS times than patients with the lowest ERBB1/EGFR expression level. High expression levels of EGFR/ERBB1 were associated with significantly increased hazard ratios for unfavorable PFS and OS outcomes in both univariate and multivariate Cox proportional hazards models. The impact of high EGFR/ERBB1 expression on the PFS and OS outcomes remained significant even after accounting for the prognostic effects of other covariates. These results regarding the prognostic effect of EGFR/ERBB1 expression were validated using the MMRF-CoMMpass RNAseq dataset generated in patients treated with more recently applied drug combinations included in contemporary induction regimens. Our findings provide new insights regarding the molecular mechanism and potential clinical significance of upregulated EGFR/ERBB1 expression in MM.


Introduction
Multiple myeloma (MM) is the second most common B-lineage lymphoid malignancy [1,2]. Recent therapeutic advances, including the incorporation of immunomodulatory drugs (IMIDs) and proteasome inhibitors (PI) into frontline and salvage regimens, have significantly improved the survival outcome of multiple myeloma (MM) [2], but treatment failures and relapse are still common due to inherent and/or acquired cancer drug resistance coupled with an immunosuppressive tumor microenvironment facilitating the immune evasion by drug-resistant MM clones [2][3][4][5][6][7][8][9][10][11][12]. MM patients with triple-class-refractory disease whose malignant clones are characterized by resistance to PIs, IMiDs, and monoclonal antibodies have an overall survival of <6 months emphasizing the urgency of identifying effective strategies to overcome triple-class cancer drug resistance [3,4]. Due to the disappointingly poor outcome of MM patients following a relapse on contemporary frontline regimens, effective treatment strategies for relapsed/refractory MM are urgently needed [1,2,4,5,8]. Personalized therapy platforms, such as precision medicines, bispecific antibodies, CAR-T cells, and antibody therapeutics continue to contribute to improved survival outcomes and thereby change the therapeutic landscape for MM [2,[6][7][8][9]. Promising clinical data with anti-CD38 monoclonal antibodies daratumumab and isatuximab and the anti-signaling lymphocyte activation marker F7 (SLAMF7) antibody elotuzumab, CD3engaging bispecific antibodies redirecting T-cells to MM targets such as CD38, the orphan G protein-coupled receptor GPRC5D, or the B-cell maturation antigen (BCMA)/CD269, as well as chimeric antigen receptor (CAR)-T cells or natural killer (NK) cells have contributed to a renewed hope regarding the potential of overcoming cancer drug resistance in MM [2]. Nevertheless, these innovative platforms have their own shortcomings causing resistance of MM clones and/or potentially life-threatening side effects, such as cytokine release syndrome (CRS) [2].
Several putative driver mutations in MM cells are considered druggable and some have been targeted with serine/threonine kinase inhibitors developed as precision medicines [2]. Examples include the development of the CDK4/6 inhibitor Palbociclib for MM with CCND1 and CDKN2C mutations [13], the development of the RAF kinase inhibitor Encorafenib (LGX818; RAF kinase inhibitor) and the MEK inhibitor Binimetinib (MEK162) for MM with BRAF V600E/K mutation [2,13,14]. However, there is insufficient information regarding the clinical potential of tyrosine kinase inhibitors (TKI) as possible components of future personalized treatment strategies for cancer drug-resistant MM.
Protein tyrosine kinases (PTK) play critical roles in normal lymphohematopoiesis, and they have also been implicated as oncoproteins in the development of B-lineage lymphoid malignancies, including leukemias, lymphomas, and MM [15][16][17]. The main goal of the present study was to identify abundantly expressed tyrosine kinases of MM cells as potential therapeutic targets with an emphasis on the relative gene expression levels of for 21 PTK, including ERBB1/epidermal growth factor receptor (EGFR), ERBB2, ERBB3, JAK1, JAK2, JAK3, TYK2, FGR, FLT3, FYN, HCK, LCK, LYN, MERTK, SRC, BLK, BMX, BTK, PTK2, SYK, TEC. The main finding of our study was differentially upregulated expression of EGFR/ERBB1 in MM cells vs. normal plasma cells which was associated with upregulated and correlated expression of several DNA binding proteins, which provides a cogent explanation for the observed EGFR/ERBB1 upregulation. The previously unknown differentially augmented expression of EGFR/ERBB1 in MM cells, especially in highrisk patients, suggests that it could serve as a therapeutic target for already approved EGFR/ERBB1 inhibitors. Intriguingly, transcript-level high EGFR/ERBB1 expression was associated with unfavorable PFS and OS outcomes in both univariate and multivariate Cox proportional hazards models. Therefore, high ERBB1/EGFR expression may also have clinical potential as a prognostic biomarker.

Data Normalization for Multiple Myeloma (MM) Samples
We downloaded raw CEL files from 3 datasets deposited in the NCBI repository (GSE171739, GSE19784, and GSE26760 from https://www.ncbi.nlm.nih.gov/geo/ (accessed on 24 September 2022)) to create a working database. We performed a probeset level normalization procedure to enable comparisons of the expression levels of specific genes in CD138 + purified malignant plasma cells from high-risk MM patients versus normal CD138 + purified normal plasma cells, as well as CD138 + , purified malignant plasma cells from standard-risk MM patients. The risk categories were based on the patients' EMC-92/ SKY-92 gene signature profiles [18][19][20][21]; we used the published SKY-92 scores to identify high-risk and standard-risk MM patients [18] in our working database. Our working database included normalized data on purified normal plasma cell samples from healthy volunteers (N = 7; GSE171739), purified malignant plasma cell samples from newly diagnosed high-risk MM patients (N = 63; GSE19784) and purified malignant plasma cell samples from newly diagnosed standard-risk MM patients (N = 219; GSE19784). GSE26760 is comprised of samples obtained from the Multiple Myeloma Research Consortium (MMRC) reference collection that were used solely for initial data normalization and quality control of the working data base.
The raw CEL files were pre-processed for batch normalization by utilizing Aroma Affymetrix statistical packages (authored and maintained by Henrik Bengtsson, Dept of Epidemiology & Biostatistics, UCSF (prev. Dept of Statistics, UC Berkeley); aroma.affymetrix_3. 2.0, aroma.core_3.2.2, aroma.light_3.24.0, affxparser_1.66.0) ran in R-studio environment (RStudio 2021.09. running with R version 4.1.2 (R Foundation for Statistical Computing, Vienna, Austria. (1 November 2021)) as previously described [22]. The PM signals were quantified using Robust Multiarray Analysis (RMA) in a 3-step process including RMA background correction, quantile normalization, and summarization by a log additive model of probes in a probeset across these samples (RmaPlm method adapted in Aroma Affymetrix). All expression values were log 2 -scaled for visualization in cluster figures. Sample annotations and patient group assignments were obtained from data files in the GEO repository (https://www.ncbi.nlm.nih.gov/geo/, accessed on 24 September 2021): GSE19784_series_matrix.txt.gz, GSE26760_series_matrix.txt.gz and GSE171739_series_matrix. txt.gz. and then accessed using the programming utility, GEOquery_2.62.1 and stringr_1.4.0, implemented in the R environment for samples with SKY-92 score determinations. Gene symbol annotation for each probeset was obtained from the database provided by the Bioconductor repository for R software (http:// www.bioconductor.org/ (accessed on 8 July 2022)) (hgu133plus2.db).
Expression values calculated from signals on the Affymetrix platform (log 2 RMA) were compared to the ranked mean level of the 100 least expressed transcripts across the 7 control samples. The mean and standard deviation values were used to determine the upper 95% confidence level to define the detection limit for the presence of specific transcripts in the MM samples. Box-plots superimposed on the dot plots were used to visualize the gene expression levels in the MM samples compared to the defined detection limit for the two EGFR/ERBB1 probesets: EGFR/ERBB1_1565483_at and EGFR/ERBB1_1565484_x_at.
Mixed Model ANOVAs implemented in the R version 4.1.2 (1 November 2021) ran in RStudio 2021.09.0 Build 351 environment, as previously described [22,23]. We investigated differential gene expression levels in purified CD138 + malignant plasma cells high-risk and standard-risk MM patients as well as purified normal plasma cells from healthy subjects. The statistical model (built using lmerTest_3.1-3 and lme4_1.1-27.1 statistical packages implemented in RStudio) controlled for variation from 2 fixed factors for "probeset" and "diagnostic group". The variation from the interaction term, "probeset × diagnostic group", was used to calculate the least square means and standard error values to pairwise comparisons of malignant vs. normal plasma cells and malignant plasma cells from high-risk vs. standard-risk MM patients. Gene chip-togene chip variation was accounted for by the random factor optimized utilizing the REML criterion (fits a variance-component model by residual (or restricted) maximum likelihood). The error term calculated for the interaction term was utilized to calculate statistical significance for linear contrasts performed between comparison groups for each probeset using emmeans_1.7.0 and lsmeans_2.30-0 statistical packages. For each set of comparisons for 65 probesets, we calculated the FDR and the 0.05 significance level and then adjusted the significance level to control for FDR [24]. p-values less than 0.05 and FDR less than 0.1 were deemed significant.
To determine the signal detection threshold for Affymetrix probesets, expression values calculated from signals on the Affymetrix platform (log 2 RMA) were compared to the ranked mean level of the 100 least expressed genes across the 7 control samples. The mean and standard deviation values were used to determine the upper 95% confidence level to define the detection limit for the expression of genes in the MM samples. Box-plots superimposed on the dot plots were used to visualize the gene expression levels in the MM samples compared to the defined detection limit for gene expression using the two EGFR/ERBB1 probesets: EGFR/ERBB1_1565483_at and EGFR/ERBB1_1565484_x_at.

Hierarchical Clustering Analysis
We used a two-way hierarchical clustering technique to organize expression patterns such that sample and probesets displaying similar expression profiles were grouped together using the average distance metric (default Euclidean distance implemented using the heatmap.2 function in the R package gplots_3.1.1), as previously described in detail [25,26].

Overall Survival and Progression-Free Survival Analysis
Archived survival data from newly diagnosed, transplant-eligible patients with newly diagnosed MM from the Dutch-Belgian Cooperative Trial Group for Hematology Oncology Group-65/German-speaking Myeloma Multicenter Group-HD4 (HOVON-65/GMMG-HD4) randomized clinical trial (ISRCTN64455289) were combined with gene expression data (GSE19784) to evaluate the impact of ERBB1/EGFR gene upregulation on progression-free survival (PFS) and overall survival (OS) of the MM patients. Patients enrolled in the HOVON-65/GMMG-HD4 study received vincristine-based (viz.: 3 cycles of vincristine 0.4 mg intravenously on days 1-4, doxorubicin 9 mg/m 2 intravenously on days 1-4, and dexamethasone 40 mg orally on days 1-4, 9-12, and 17-20 or bortezomib-based (viz.: bortezomib 1.3 mg/m 2 was given intravenously on days 1, 4, 8, and 11, doxorubicin 9 mg/m 2 intravenously on days 1-4, and dexamethasone 40 mg orally on days 1-4, 9-12, and 17-20) standard induction chemotherapy After induction therapy, patients received one  or two (GMMG-HD4) cycles of high-dose melphalan (200 mg/m 2 intravenously) with autologous stem-cell rescue followed by maintenance treatment with thalidomide (50 mg per day orally; group assigned to vincristine-based induction treatment) or bortezomib (1.3 mg/m 2 intravenously once every 2 weeks; group assigned to bortezomibbased induction treatment) for 2 years. RMA-normalized values from 282 newly diagnosed MM patients (pooled standard-risk and high-risk samples from GSE19784) were rank ordered according to the expression of the ERBB1_1565484_x_at probeset. PFS and OS times were evaluated for 280 MM patients with available follow-up data (data obtained from the supplementary file 41375_2012_BFleu2012127_MOESM30_ESM.xls in from [18] that enabled comparisons of OS and PFS for these MM patients with the highest expression of ERBB1/EGFR (top 40th percentile; N = 112) versus MM patients with the lowest expression of ERBB1/EGFR (bottom 40th percentile; N = 112). Out of the 112 MM patients that expressed the highest level of ERBB1/EGFR, 33 were categorized as high-risk using the SKY-92 classifier (29%), in comparison to the 112 MM patients categorized as 18 high-risk (16%) expressing the lowest level of ERBB1/EGFR. The Kaplan-Meier (KM) method was used to investigate the PFS and OS outcome data utilizing software packages survival_3.2-13, survminer_0.4.9 and survMisc_0.5.5 operated in the R environment. Statistical significance between comparison groups was determined using the Log-rank chi-squared values and p-values less than 0.05 were deemed significant.
Multivariate analysis of the effect of EGFR/ERBB1 expression on the PFS and OS outcomes was performed using the Cox proportional hazards model to adjust for patient prognostic characteristics as variables that included the prognostic staging data according to the International Staging System (ISS), and the 92-gene SKY signature assessed risk score. We tested whether the correlation of expression level (log 2 RMA) was a significant independent prognostic factor controlling for each as well as all of the other variables. Estimates of the life table hazard ratios (HR) were calculated using the exponentiated regression coefficient for proportional hazards model analyses implemented in R (survival_3.2-13 ran in R version 4.1.2 (1 November 2021)). Forest plots were utilized to visualize the HRs for Cox proportional hazards model (survminer_0.4.9 ran in R version 4.1.2 (1 November 2021)).

Processing of the Multiple Myeloma Research Foundation (MMRF)-CoMMpass RNAseq Dataset
The entire database from 33 TCGA projects has been migrated to the NCI's Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov/, accessed on 1 October 2022). One of the TCGA projects, the Multiple Myeloma Research Foundation (MMRF) CoMMpass study, harbors clinical information for 995 cases and RNAseq data files for 787 cases and are all made available via the GDC portal. We used the MMRF-CoMMpass RNAseq dataset for validation of the microarray-based findings regarding the prognostic effect of EGFR/ERBB1 expression. The workflow for processing RNAseq datasets standardizes the sequence alignment to the same genome build (GRCh38.p0) and harmonized quantification techniques (STAR-count; detailed in https://docs.gdc.cancer.gov/ Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/, accessed on 20 September 2022) to enable meta-analysis across multiple projects.
We utilized the latest version of the RNAseq STAR-quantified data for all 787 cases (release date 29 March 2022, version 32, STAR is implemented as a standalone C++ code. STAR is free open source software distributed under GPLv3 license and can be downloaded from http://code.google.com/p/rna-star/, accessed on 20 September 2022. And analysis performed by GDC) whereby the STAR-aligned read groups using the two-pass method were employed to generate the final alignments to calculate the expression metrics, and the quality of sequence reads was assessed by best matching the alignments of the reads to the reference genome sequence using pre-alignment with FASTQC and post-alignment with Picard Tools.
Gene level RNAseq raw count data (STAR-aligned unstranded number of reads aligned per gene per sample) were downloaded from the archived MMRF CoMMpass dataset by connecting to the GDC portal: https://gdc.cancer.gov/about-gdc/contributed-genomicdata-cancer-research/foundation-medicine/multiple-myeloma-research-foundation-mmrf. (accessed on 24 September 2022) using the Bioconductor package, Genomic-DataCommons_1.18.0 (Release date 2021-08-11 with full functionality as provided by TCGAbiolinks for accessing GDC data; https://bioconductor.org/packages/release/bioc/ html/GenomicDataCommons.html (accessed on 20 September 2022)) implemented in R version 4.1.2 (1 November 2021). The mRNA expression data were deposited in files appended with " . . . .rna_seq.augmented_star_gene_counts.tsv". Clinical data for each MM patient were also acquired by functions provided in GenomicDataCommons_1.18.0 and the case IDs were matched with RNAseq unique identifiers utilizing the metadata ("metadata.cart.2022-09-04.json") from the GDC portal converted into R data by running the utilities rjson_0.  [27]. For each gene, DESeq2 uses a generalized linear model that fits read counts to negative binomial distribution to calculate the mean and variance estimates, whereby the mean is taken as a quantity proportional to the concentration of cDNA fragments from the gene in the sample and scaled by a normalization factor across all samples. The normalization method in the DESeq2 algorithm determines the counts divided by sample-specific size factors calculated from the median ratio of gene counts relative to the geometric mean per gene across all samples that accounts for sequencing depth and RNA composition for each gene. This method allows for fold-change comparisons across treatment groups in the GLM model [28][29][30]. Statistical significance was assessed by testing the null hypothesis that there is no differential expression across the two sample groups (Log 2 fold change = 0) using the Wald test [27] reporting the test statistic and p-value for each gene. Genes were considered differentially expressed if the p-values were less than 0.05 after adjusting for multiple comparisons [24]. To visualize the gene expression levels in boxplots, the normalized log 2 values were calculated from the RNAseq count data using the variance stabilization method supplied by the algorithms in vsn_3.62.0 [31]. Low count values tend to generate large fold-changes. Therefore, to calculate a more accurate log 2 fold-change estimate, we applied a shrinkage of the log fold-change estimates toward zero when the read counts were low and variable ("normal" function in the DESeq2 package") [32].

Overall Survival Analysis and Cox Proportional Hazards Model of the MMRF CoMMpass Study
The Kaplan-Meier (KM) method, log-rank chi-square test, was used to investigate the PFS and OS outcomes of MM patients utilizing the software packages survival_3.2-13, survminer_0.4.9 and survMisc_0.5.5 operated in the R environment. Statistical significance between comparison groups was determined using the log-rank chi-squared values and p-values less than 0.05 were deemed significant. Graphical representations of the survival curves were visualized using graph-drawing packages implemented in the R programming environment: dplyr_1.0.7, ggplot2_3.3.5 and ggthemes_4.2.4. OS curves were compared for patients expressing high levels of EGFR/ERBB1 (top 50th percentile) versus low levels of EGFR/ERBB1 (bottom 50th percentile). PFS data and biochemical measurements (beta 2 microglobulin and albumin levels) were obtained from the web repository, doi.org/10.6084/ m9.figshare.11941494 [33]. The percentiles of patients expressing EGFR/ERBB1 were determined using the metric fragments per kilobase of transcript per million mapped reads (FPKM) calculation normalized to the upper quartile FPKM (FPKM-UQ) for the whole geneset in each sample.
Multivariate analysis of the effect of MMRF CoMMpass cohort groups on the PFS and OS outcomes was performed using the multivariate Cox proportional hazards model to adjust for patient prognostic characteristics. The model included (i) the expression level (high versus low FPKM-UQ values) of the EGFR/ERBB1 gene, (ii) the prognostic staging data according to the International Staging System (ISS), (iii) age, (iv) gender, (v) serum Beta 2 microglobulin levels and (vi) serum albumin levels. The multivariate Cox proportional hazards model was compared to univariate models for each of the variables to assess the independent effect of high levels of EGFR/ERBB1 expression on the PFS and OS outcomes in MM patients. Estimates of the life table hazard ratios (HRs) were calculated using the exponentiated regression coefficient for Cox proportional hazards analyses implemented in R (survival_3.2-13 ran in R version 4.1.2. Forest Plots were utilized to visualize the Hazard ratios for Cox proportional hazards model (survminer_0.4.9 ran in R version 4.1.2 (1 November 2021)).

Amplified Expression of Selective Tyrosine Kinase Genes in Malignant Plasma Cells from MM Patients
We first compared the transcriptomes of purified CD138 + malignant plasma cells from 282 newly diagnosed MM patients with the transcriptomes of CD138 + normal plasma cells from 7 healthy volunteers. Epidermal growth factor receptor (EGFR) ERB1, ERBB3, SRC, and MERTK were significantly overexpressed in MM cells in comparison to normal control plasma cells (Figure 1, Table S1). ERBB1_1565483_at was the most significantly upregulated probeset (Fold Change = 16.86; p-value < 1 × 10 −8 ) followed by ERBB1_1565484_x_at (Fold Change = 9.91; p-value < 1 × 10 −8 ) and ERBB1_211607_x_at (Fold Change = 2.54; p-value = 2.7 × 10 −6 ) (Table S1). It is noteworthy that the absolute EGFR/ERRB1 gene expression levels from all the MM samples regardless of the risk category or ISS stage were above the lowest expressed transcripts on the Affymetrix Human Genome U133 Plus 2.0 platform for both probesets ( Figures S1 and S2, Tables S2 and S3).
Five probesets for ERBB1/EGFR and 2 probesets for ERBB3 were significantly overexpressed in MM cells. SRC (2 probesets), ERBB1/EGFR, ERBB3 (2 probesets) and MERTK formed a cassette of upregulated probesets in these pooled MM patients. The gene for the MM surface antigen BCMA, selected as a control marker gene for MM, was expressed at a 2.4-fold higher level in malignant plasma cells from MM patients than normal plasma cells (Figure 1, Table S1). BCMA was co-regulated with 2 probesets for ERBB1/ EGFR.

Differentially Amplified Expression of ERBB1/EGFR and JAK3 Genes in Malignant Plasma Cells from High-Risk MM Patients
We compared the transcript-level PTK expression profiles of malignant plasma cells from high-risk (N = 63) vs. standard-risk (N = 219) MM patients. Notably, the expression levels of ERBB1/EGFR and JAK3 were differentially amplified in malignant plasma cells from high-risk MM patients ( Figure 3, Table S5).

Amplified Expression of ERBB1/EGFR Is Associated with Poor PFS and OS in Newly Diagnosed MM
We next examined the effects of ERBB1/EGFR expression on PFS and OS in univariate and multivariate Cox proportional hazards models that also included SKY-92 gene signaturebased risk assignment and prognostic ISS stage as categorical covariates. High EGFR/ ERBB1 expression was associated with significantly increased hazard ratios (HR) for unfavorable PFS outcomes in both univariate (HR = 1.85 (1.3-2.62)) and multivariate (HR = 1.5 (1.05-2.15)) ( Figure S5A,B). Cox proportional hazards models. High EGFR/ ERBB1 expression was also associated with significantly increased HRs for unfavorable OS outcomes in both univariate (HR = 2.28 (1.38-3.78)) and multivariate (HR =1.75 (1.05-2.92)) Cox proportional hazards models ( Figure S5C,D). The impact of high EGFR/ERBB1 expression on the PFS and OS outcomes remained significant even after accounting for the prognostic effects of the risk category and ISS stage (Supplementary Figure S5A-D).

Analyses Using the MMRF-CoMMpass RNAseq Dataset for Validation of the Microarray-Based Findings from the HOVON65/GMMG-HD4 Trial Regarding the Prognostic Effect of EGFR/ERBB1 Expression
The randomized HOVON65/GMMG-HD4 Phase 3 study evaluating the clinical efficacy of bortezomib induction/maintenance versus classic cytotoxic drugs and thalidomide maintenance relied in part on cytotoxic agents such as doxorubicin and vincristine that are no longer routinely used as part of frontline induction therapy in MM patients [50].
To ensure that EGFR/ERBB1 expression remains relevant with modern induction regimens, we next used the MMRF-CoMMpass RNAseq dataset generated in patients treated with contemporary backbone combinations of immunomodulatory drugs (IMiDs), proteosome inhibitors, and dexamethasone for validation of the microarray-based findings from the HOVON65/GMMG-HD4 study regarding the prognostic effect of EGFR/ERBB1 expression [51][52][53]. We first determined the lowest detection level for the normalized, variance stabilized log 2 count value at which there were zero alignments to the EGFR/ERBB1 gene in 766 MM patients with ISS staging information. Only cells from 157 patients (20.5%) across all 3 prognostic ISS stages had zero alignments to the EGFR/ERBB1 gene ( Figure S6).
The Kaplan-Meier (KM) method, log-rank chi-square test, was then used to investigate the PFS (N = 598) and OS (N = 716) of evaluable MM patients in relationship to their RNAseq-based EGFR/ERBB1 expression levels. Patients with higher levels of EGFR/ ERBB1 expression had worse outcomes with significantly shorter PFS and OS times ( Figure  S7). We next examined the effects of ERBB1/EGFR expression on PFS and OS in univariate and multivariate Cox proportional hazards models. Higher EGFR/ERBB1 expression was associated with significantly increased hazard ratios (HR) for unfavorable PFS outcomes in both univariate (HR = 1.85 [1.3-2.62]) and multivariate (HR = 1.5 [1.05-2.15]) ( Figure  5A) Cox proportional hazards models. Higher EGFR/ERBB1 expression was also associated with significantly increased HRs for unfavorable OS outcomes in both univariate (HR = 2.28 (1.38-3.78)) and multivariate (HR =1. 75 (1.05-2.92)) Cox proportional hazards models ( Figure 5B). The impact of higher EGFR/ERBB1 expression on the PFS (HR: 1.37, p = 0.005) and OS (HR: 1.73, p = 0.017) outcomes remained statistically significant even after accounting for the prognostic effects of the prognostic ISS stage, age, gender, serum albumin level and serum beta 2 microglobulin level ( Figure 5).
Gene level normalized expression values (FPKM-UQ) obtained from GDC data portal aligned to the EGFR/ERBB1, sequences (GRCh38.p0 genome build)) was correlated with PFS ( Figure 5 Panels A, B) and OS ( Figure 5 Panels C, D) times using the univariate and multivariate Cox proportional hazards models. We investigated 3

Discussion
Our detailed comparative analysis of the transcript-level PTK expression profiles of purified CD138 + malignant plasma cells from newly diagnosed MM patients vs. normal plasma cells showed that ERBB1/EGFR, ERBB3, SRC, and MERTK were significantly overexpressed in MM cells in comparison to normal control plasma cells. Further, ERBB1/EGFR and JAK3 exhibited differentially amplified expression in malignant plasma cells from high-risk MM patients. These results provide new insights regarding the clinical impact potential of targeting ERBB1/EGFR in high-risk MM.
Notably, amplified expression of ERBB1/EGFR emerged as a prognostic biomarker in both the dataset from the HOVON65/GMMG-HD4 study (depicted in Figure S4) as well as the MMRF-CoMMpass RNAseq validation dataset (depicted in Figure S7): Patients with the highest ERBB1/EGFR expression levels had significantly worse PFS and OS outcomes than patients with lowest ERBB1/EGFR expression levels. Hence, this finding was not an artifact unique to a specific clinical study. Furthermore, the impact of high EGFR/ERBB1 expression on the PFS and OS outcomes remained statistically significant even after accounting for the effects of the SKY-92 gene signature-based risk assignment, prognostic ISS stage, age, gender, serum albumin level and serum beta 2 microglobulin level in multivariate Cox proportional hazards models shown in Figures S5 and 5.
Tumors with activating ERBB1/EGFR exon 19 deletions and exon 21 mutations have exhibited favorable responses to TKI of ERBB1/EGFR [58]. Although the first generation TKI of ERBB1/EGFR have not shown meaningful clinical activity in cancer patients with wildtype ERBB1/EGFR or mutant ERBB1/EGFR with T790M, new-generation TKI such as Afatinib are active even against wildtype or T790M mutant ERBB1/EGFR [56]. It is noteworthy that tumors with ERBB1/EGFR exon 20 insertion mutations are Uckun and Qazi Page 13 significantly less sensitive to both TKI and monoclonal antibodies targeting ERBB1/EGFR [62]. Clinical evaluation of ERBB1/EGFR targeting therapeutics in MM patients would benefit from integrated biomarker research to determine if MM cells display or acquire any ERBB1/EGFR mutations that may affect treatment outcomes. TKI of ERBB1/EGFR cause treatment-emergent side effects some of can be associated with potentially fatal serious toxicities such as lung injury/interstitial lung disease, especially in patients with preexisting pulmonary disease [63,64]. Careful monitoring and risk mitigation strategies will be required to maximize patient safety in clinical trials of ERBB1/EGFR inhibitors in MM patients.
Our results extend previous work that has implicated the EGF-family growth factor AREG as an autocrine growth factor that binds selectively to ERBB1/EGFR on MM cells and promotes their growth and survival [65]. Further, the ERBB1/EGFR ligand EGF derived from bone marrow stromal cells has been shown to stimulate the clonogenic growth of MM cells [66]. The ERBB1/EGFR inhibitor Iressa induced apoptosis in primary MM cells and enhanced the apoptotic activity of dexamethasone [65]. RAS or RAF mutations may render MM cells resistant to treatment with ERBB1/EGFR1 inhibitors and limit any potential clinical benefit to MM patients with a KRAS/NRAS/BRAF triple wildtype genotype [67]. MM-derived exosomes have been shown to contain AREG and activate the EGF signaling pathway in the bone microenvironment leading to osteoclastogenesis [68].
As AREG selectively binds to ERBB1/EGFR, TKI of ERBB1/EGFR could also inhibit the progression of osteolytic bone disease in MM. It is noteworthy that the anti-EGFR antibody Cetuximab was evaluated in a pilot study as a single agent in patients with refractory or relapsed MM who had previously received at least one line of prior treatment and were not eligible to undergo autologous stem cell transplantation [69]. VonTreschkow et al. reported an overall response rate of 27% with acceptable tolerability [69].
We hypothesized that the augmented expression of the EGFR/ERBB1 in MM cells could be owing to the amplified expression of transcription factors that are known to upregulate EGFR/ERBB1 expression. We therefore examined if the ERBB1/EGFR expression in MM cells is correlated with expression levels of such DNA binding proteins. Our analyses provided unprecedented evidence that the amplified expression of EGFR/ERBB1 mRNA in MM cells was positively correlated with increased expression levels of mRNAs for several DNA binding proteins/transcription factors, including ETF [34][35][36], SP-1 [37][38][39][40][41], TCF [42,43], HOXB5 [44,45], RPF-1 [46] and AP-1 [33,[47][48][49]. Based on these findings we propose a model according to which the mechanism of the observed transcript-level upregulation of the EGFR/ERB1 gene in MM samples is related to the upregulated expression of the correlated DNA binding proteins ( Figure 6). The mRNA for these transcription factors were highly positively correlated with EGFR/ERBB1 gene upregulation (p < 0.0001) in MM patients.
Correlation coefficients of gene expression levels calculated between EGFR/ERBB1 and transcription factor probesets (282 MM patients (GSE19784)) identified 11 genes that were significantly positively correlated with EGFR/ERBB1 expression. Depicted in Figure 6 are the mRNAs that coded for transcription factor proteins (oval shapes) bound to DNA sequence binding sites for these transcription factors (rectangles).
Janus family protein tyrosine kinases, including JAK1, JAK2, JAK3 and TYK2, have important regulatory roles in multiple signal transduction pathways intimately linked to survival, proliferation, maturation, and function of normal as well as malignant lymphoid cells [86][87][88]. An intriguing new finding gained in the present study was the differential upregulation of JAK3 expression in malignant plasma cells from high-risk MM patients. JAK3 associates with the common cytokine receptor γ chain (γc) and thereby regulates signaling by the cytokines Interleukin (IL)-2, IL-4, IL-7, IL-9, IL-15, and IL-21 [89]. Tofacitinib (Xeljanz ® ) is a potent, selective JAK inhibitor that preferentially inhibits JAK1 and JAK3 used for the treatment of moderate to severe active rheumatoid arthritis (RA) ( Table 1) [90][91][92], Tofacitinib had promising in vitro activity on MM cells [93]. Further, the expression of JAK1 and JAK2 were previously reported to be upregulated in some MM cells and the JAK1/JAK2 inhibitor Ruxolitinib in combination with Bortezomib, Itacitinib or Daratumumab inhibited JAK/STAT3 phosphorylation, inhibited in vitro and in vivo myeloma cell growth and induced cell apoptosis [94][95][96]. In a phase I clinical trial, Ruxolitinib exhibited promising early clinical activity in R/R MM patients [97], Contemporary induction protocols for MM employ Pis such as bortezomib and Carfilzomib [2], PIs trigger the apoptotic death of MM cells by contributing to an exaggerated unfolded protein response (UPR) pathway and ER stress [2]. They further impair the survival-promoting interactions between MM cells and stromal elements in the bone marrow microenvironment, and they promote the immunogenic death of damaged MM cells [2]. The combination of PIs with JAK inhibitors may potentiate their pro-apoptotic activity against MM cells. Likewise, JAK inhibitors could be combined with other apoptosis-inducing MM drugs, such as the BCL-2 homology 3 (BH3)-mimetic Venetoclax and inhibitors of MCL-1 such as AMG-176 and MIK665 [1]. Recently, serious concerns have emerged regarding the safety profile of JAK inhibitors, especially JAK3 inhibitor Tofacitinib as well as Baricitinib,and Upadacitinib [98]. In particular, there was an increased risk of serious heart-related events such as heart attack or stroke, cancer, blood clots, and death with arthritis and ulcerative colitis medicines Xeljanz and Xeljanz XR (tofacitinib). Therefore, targeting JAK3 in MM will require risk mitigation strategies to reduce the incidence of serious cardiovascular side effects. Several new JAK3 inhibitors have been developed and some of these new-generation inhibitors could have more favorable safety profiles [99][100][101][102]. The favorable safety profile and promising anti-tumor activity of the dual-function JAK3-ERBB1/EGFR inhibitor WHI-P131 also warrants further studies [103].

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.
Five out of the 10 upregulated probesets were for the ERBB1 gene and 2 probesets for the ERBB3 gene. In addition, MERTK, SRC and BCMA were also significantly upregulated in Uckun and Qazi Page 25 malignant plasma cells. Co-regulation of BCMA and 2 probesets of ERBB1 are depicted in the cluster figure. The unfavorable impact of high EGFR/ERBB1 expression on PFS and OS outcomes of MM patients from the MMRF-CoMMpass study in univariate and multivariate Cox proportional hazards models. (A,B) Comparison of the PFS times for 582 evaluable patients showed a significant increase in HR for patients with high expres-sion of EGFR/ERBB1 compared to patients with low EGFR/ERBB1 expression in the multivariate model. (C,D) OS relationships were evaluable for 695 patients). Significant p-values in B and D are indicated by *, ** and *** for p < 0.05, p < 0.01 and p < 0.001 respectively. Transcription factor genes that are positively correlated with EGFR/ERRB1 expression in Multiple Myeloma (MM) patients.