The microarray mRNA expression profiles and related clinical information for GSVA analysis were obtained from GSE24080 on the GEO database (https://www.ncbi.nlm.nih.gov/geo/). External validation GEO datasets include GSE7982 and GSE2658. The transcriptome of different developmental stages of MM was obtained from GSE6477. Transcriptome data of sPCL samples were obtained from GSE39683 and GSE66291. The TCGA-MMRF-COMPASS project (https://portal.gdc.cancer.gov/projects) consisted of 787 MM transcriptome data along with corresponding clinical information. Oncomine datasets (https://www.oncomine.org/resource/login.html) include Agnelli myeloma 3, Zhan myeloma 2. The single-cell transcriptome data were obtained from GSE106218. The expression matrix of GEO, TCGA datasets along with detailed clinical characteristics of COMPASS-MMRF patients used in this study can be accessed and cited via Zheng, Bo (2021), “The distinctive role of YAP1 in multiple myeloma”, Mendeley Data, V2, doi: 10.17632/mmp7cw7rx9.2.

The GSVA software package (V1.25.4) for R was applied as a non-parametric, unsupervised method for estimating the variation of key gene sets in MM. The input for the GSVA algorithm was a gene expression matrix of log2 microarray expression values and a collection of C6-oncogenic signature gene sets from the Molecular Signature Database (MSigDB, version 7.4). GSVA scores were determined nonparametrically using a Kolmogorov–Smirnov (KS)-like random walk statistic and a negative value for a particular sample and gene set.

The “survival” and “survminer” R packages were applied for survival analysis. Surv_cutpoint function was applied to acquire the optimal cutoff value. The survival information of GEO and TCGA datasets was acquired along with their expression matrix data. p-value less than 0.05 was considered significant. A single-factor Cox model was used to determine whether a single gene was related to the prognosis of MM.

The “Clusterprofiler” package (V3.16.0) of R language was applied for enrichment analysis of Reactome pathways. “Clusterprofiler” is an R package of Bioconductor, which can perform statistical analysis and visualization of functional clustering on gene sets or gene clusters. When the adjusted p-value was below 0.05, the Reactome pathways were identified as significantly enriched by these genes.

To unveil biological correlations of the obtained gene expression profiles, the transcriptome data were compared using GSEA (http://www.gsea-msigdb.org/gsea/). GSEA uses a weighted Kolmogorov–Smirnov method to determine whether the distribution of genes in the gene set is different from the normal distribution. h.all.v7.4.symbols and c2.cp.reactome.v7.4.symbols were selected for analysis. False discovery rate (FDR) < 0.05 and adjusted p-value < 0.05 were considered statistically significant.

“Monocle2 (v2.16.0)”, an R package, was applied to conduct single-cell evolutionary trajectory analysis and to estimate the transcriptome evolution in primary and extramedullary MM samples. Monocle is based on the assumption that one-dimensional “time” can depict the multi-dimensional expression values to elucidate the cell state transitions. In the trajectory analysis, we used genes meeting the following standards: mean_expression ≥ 0.1 and dispersion_empirical ≥ 1 * dispersion_fit to sort cells in pseudo-time order. The visualization functions “plot_cell_trajectory” were used to plot the minimum spanning tree on cells. Genes that changed along with the pseudotime were calculated (q-value <0.01) by the “differentialGeneTest” function and visualized with the plot_pseudotime_heatmap and the genes were clustered into subgroups according to the gene expression patterns.

All statistical analyses were performed using R 4.1.0 and GraphPad prism 8. Two-sided Student’s t-test for unpaired samples was applied to evaluate the significance of differences in experiments. Pearson correlation coefficients were calculated to measure associations among the mRNA expression level of various genes. OS and RFS were assessed with the Kaplan–Meier method, and the differences between the groups were compared by the log-rank test. Values of p and q less than 0.05 were considered statistically significant. Other key packages used in this study include GGally, ggrisk.

We performed GSVA analysis on the GEO dataset (GSE24080), which contains the transcriptome data of bone marrow purified plasma cells from 559 newly diagnosed MM patients. The median age of this cohort is 58 and male accounts for 60.3% ( Table S1 ). GSVA analysis aims to investigate whether a set of genes are randomly distributed in a specific phenotype compared with others and identifies biological or pathological correlations. To this end, we examined the expression of C6-oncogenic signature gene sets from the Molecular Signature Database (MSigDB, version 7.4) in all 559 MM cases. The result showed that 599 patients were separated into two groups ( Figure S1A ). We further identified 5 signaling pathways that fit the following criteria: significantly up/downregulated in one MM Group while exhibiting the opposite in the other group. In the heatmap ( Figure 1A and Table S2 ), we observed that in Group 1, which contains 231 patients, VEGFA and ERBB2 pathways were activated while other oncogenic pathways such as MYC, mTOR, and YAP1 exhibited the opposite. This unusual phenomenon might be attributed to the different molecular landscapes between solid and blood cancers. Notably, patients in Group 1 suffered from a significantly shorter survival (p < 0.0001) and higher recurrence rate (p < 0.001) compared with Group 2 ( Figures 1B, C ). Next, we investigated the expression level of the above 5 key genes in the pan-cancer scenario in TIMER database and found that YAP1 was the only gene significantly downregulated in blood cancers while consistently highly expressed in solid tumors ( Figure S1B ).

Considering that the above 4 oncogenic pathways were specifically dysregulated along with YAP1, we moved on to investigate the correlation among them. The result showed that the YAP1 pathway was positively correlated with MYC and mTOR pathways while negatively correlated with VEGFA and ERBB2 pathways ( Figure 2A ). The Hippo/YAP1/TAZ pathway is well known as being highly conservative and the activation of the YAP1/TAZ complex is responsible for cell proliferation and anti-apoptosis (7). In solid tumors, the YAP1/TAZ complex improves the transcription of a series of oncogenes including MYC (8). In our study, YAP1 and TAZ were downregulated in Group 1 while MYC along with VEGFA and mTOR exhibited the opposite ( Figure 2B ). YAP1 showed a weak while statistically significant negative correlation with MYC and VEGFA ( Figure S2A ). Furthermore, YAP1 exhibited a positive prognostic factor in MM while MYC was a negative prognostic factor ( Figure 2C ). YAP1+MYC- patients showed the most favorable prognosis while YAP1-MYC+ patients showed the worst prognosis ( Figure 2D ). TAZ, VEGAF, and ERBB2 were all identified as negative prognostic factors in MM patients ( Figure S2B ).

The development of multiple myeloma goes through monoclonal gammopathy of undetermined significance (MGUS), smoldering multiple myeloma (sMM), symptomatic MM, and relapsed MM. Based on the data from GSE6477 and GSE5900 (9), we discovered that the expression of YAP1 decreased gradually from NC to sMM until a quick drop is experienced when progressing into MM. This phenomenon showed that YAP1 level was negatively correlated with tumor burden in MM. Comparably, MYC exhibited a gradual increase from NC to relapsed MM cases ( Figure 3A ). The same trend was validated in an external Oncomine dataset (Angelli Myeloma 3) with an exception of normal plasma cells probably due to fewer cases enrolled ( Figure S3B ). The prognostic significances of YAP1/MYC were further confirmed in 4 external independent datasets including GEO datasets (GSE9782 and GSE2658), TCGA dataset (MMRF-COMPASS), and Oncomine dataset (Zhan myeloma 2) in over 2,000 MM patients ( Figures 3B, C and Figures S3C, D ). Treatment information of 787 MM patients was acquired from the MMRF-COMPASS study ( Table S3 ). Notably, the expression level of YAP1 efficiently affected the response of various treatments. In MM patients with low YAP1 expression, the carfilzomib-based treatment showed obvious superiority over other treatments compared with other MM patients ( Figure 3D ). In MM patients with high MYC expression, the efficacy of IMiD treatment was improved compared with MM patients with low MYC expression ( Figure S3E ) and this trend was further amplified when we compared YAP-MYC+ MM patients with the rest ( Figure 3E and Figure S3F ).

It was reported that low YAP1 level in hematological cancers resulted in an improved DNA repairing process, which facilitated cell proliferation and improved oncogenicity (9). We observed that DNA repair-related genes such as PARP1 and H2AX were improved in Group 1 while apoptosis-related gene TP73 was suppressed ( Figure 4A ). According to previous studies, in extramedullary myeloma cells, the glycolysis process and oxidative phosphorylation process were enhanced (10). We found that glycolysis-related genes and oxidative phosphorylation-related genes were significantly upregulated in Group 1 patients. Based on the STRING database, the protein–protein interaction (PPI) network of the genes in Figure 4A was drawn ( Figure 4B ). The GO analysis result showed that multiple key energy pathways were involved ( Figure 4C ), indicating a metabolic reprogramming in Group 1 patients, which might lead to its increased oncogenicity and unpromising prognosis. In the expression correlation heatmap, we observed that in Group 1 patients, metabolic genes were mostly co-expressed ( Figure 4D ). The univariate analysis result showed that most metabolic genes posed as oncogenes while YAP1 and TP73 posed as tumor suppressor genes ( Figure 4E ). In GSE39683 and GSE66291 datasets, we acquired the transcriptome data of MM and second plasma cell leukemia (sPCL) and found that YAP1 was also downregulated in sPCL while MYC showed no difference ( Figure S4A ). With further exploration, we discovered that the activity of DNA repair, glycolysis, and oxidative phosphorylation was enhanced in sPCL, which simulated the phenotype of YAP1-MYC+ MM ( Figure S4B ), indicating that lower YAP1 expression might facilitate extramedullary invasion.

To further analyze the intratumoral heterogeneity in terms of YAP1/MYC expression, we collected the single-cell transcriptome data of 477 myeloma cells from GSE106218 which contained 9 MM samples and 3 paired EM samples. However, YAP1 data were not available on the single-cell level due to its low expression. Hence, we focus on the role of MYC in the process of extramedullary invasion. First, we performed dimension reduction analysis and displayed the sc-seq data on tSNE to investigate the similarities and divergence among different clusters and samples ( Figure 5A ). The result showed that 3 EM samples were well isolated from the clusters of primary MM samples except for MM17. We found that DNA repair and metabolic genes were marked upregulated in EM cases along with MM17 ( Figure 5B ). To examine the evolutionary transcriptomic change between primary MM and extramedullary myeloma, we performed a trajectory analysis for all 12 samples. Pseudotime analysis indicated a branched evolution from primary myelomas to extramedullary myelomas ( Figure 5C and Figure S5A ). Intriguingly, along the two extramedullary branches, MM02EM and MM36EM stayed together while MM34EM was located on the other branch. Comparably, MM34EM expressed the highest level of DNA repair and metabolic genes along with MYC ( Figure 5B ). Also, the MM34EM case experienced the shortest extramedullary invasion time (3 months) compared with MM02EM and MM36EM (both 20 months) ( Table S4 ). MM17, which showed a similar phenotype as EM samples ( Figure 5B ), is located on the same extramedullary branch with MM34EM, indicating a propensity for EM invasion ( Figure 5C ). Notably, MM17 had the shortest survival time (3 months) among all primary MM cases. The expressions of DNA repair and metabolic genes along with MYC were marked enhanced on the extramedullary branches and experienced a downfall from EM to MM cases ( Figures 5D, E )

In the results above, we identified a MM subgroup with a unique molecular feature, prognostic significance, and therapeutic response. Hence, we next intended to explore deeper into the difference of its mRNA expression pattern. To investigate the transcriptome signature of Group 1 (YAP1-MYC+) and Group 2 (YAP1+MYC-) cases, we firstly explored the distribution of all samples based on similarities in gene expression data. The PCA analysis revealed huge nonoverlapping areas between the 2 groups ( Figure 6A ). Gene expression analysis identified a total of 16,034 genes differentially expressed between two groups, of which 4,517 were overexpressed and 11,517 genes were underexpressed in Group 1 relative to Group 2 patients ( Figure 6B and Table S5 ). Exemplary DEGs were annotated in the volcano plot and heatmap ( Figure 6C and Figure S6A ). The list of the 50 genes with the largest and smallest fold change (FC) is shown in Figure 6D . Functional analysis using the Reactome database revealed statistically significant enrichment for mitotic prometaphase, M phase, mitotic anaphase, mitotic metaphase/anaphase, etc. ( Figure 6E and Figure S6B , Table S6 ). Since the functional analysis in Figure 6E did not take the fold change of gene expression into account, we performed GSEA analysis to identify the key changes in biological processes in terms of reactome database and Hallmark gene sets (MSigDB, version 7.4) between Group 1 and Group 2 ( Tables S7 , S8 ). The result showed that cell cycle was accelerated in Group 1 patients along with enhanced glycolysis, DNA repair, and oxidative phosphorylation processes. These findings confirmed that the improved oncogenicity in Group 1 MM cases resulted in high cellular proliferation ( Figure 6F and Figure S6C ). Additionally, we discovered that MM patients with the presence of soft tissue plasmacytoma exhibited enhanced activity of enhanced glycolysis, DNA repair, and oxidative phosphorylation processes along with mitosis ( Figure S6D ). Also, a subgroup of MMRF patients with low YAP1 expression and high MYC expression exhibited stronger activities of DNA repair/Glycolysis/Oxidative phosphorylation ( Figure S6E ).

From the above results, we concluded that DNA repair/glycolysis/oxidative phosphorylation processes were improved in Group 1, which exhibited enhanced oncogenicity and stronger propensity for extramedullary invasion. We identified that higher activities of these 3 processes were significantly correlated with worse survival in MM ( Figure 7A and Figures S7A–C ). An integrated risk score based on the above 3 processes was calculated for each MM patient using the COX regression model ( Figure 7B ). The optimal cutoff value was selected by using the Gordon index (cutoff = 0.21), and KM curves of survival was also performed in Figure 7C . The same result was confirmed in the MMRF-COMPASS dataset ( Figure S7D ). Furthermore, we performed a multivariable analysis corrected for ISS stages using the MMRF dataset. The result showed that YAP1 expression and glycolysis activity functioned as prognostic factors even in a multivariable fashion ( Figure S7E ). Also, we observed that enhanced DNA repair and glycolysis activity would impair the therapeutic response of bortezomib-based and bortezomib+carfilzomib-based therapy respectively ( Figures 7D, E ). Interestingly, we discovered that the presence of unfavorable chromosomal alteration such as t(4;14) and t(11;14) also enhanced the activity of DNA repair/glycolysis/oxidative phosphorylation processes ( Figures S7F, G ). The entire study is summarized in Figure 8 .