Prognostic significance of esterase gene expression in multiple myeloma

Background Esterase enzymes differ in substrate specificity and biological function and may display dysregulated expression in cancer. This study evaluated the biological significance of esterase expression in multiple myeloma (MM). Methods For gene expression profiling and evaluation of genomic variants in the Institute for Molecular Medicine Finland (FIMM) cohort, bone marrow aspirates were obtained from patients with newly diagnosed MM (NDMM) or relapsed/refractory MM (RRMM). CD138+ plasma cells were enriched and used for RNA sequencing and analysis, and to evaluate genomic variation. The Multiple Myeloma Research Foundation (MMRF) Relating Clinical Outcomes in MM to Personal Assessment of Genetic Profile (CoMMpass) dataset was used for validation of the findings from FIMM. Results MM patients (NDMM, n = 56; RRMM, n = 78) provided 171 bone marrow aspirates (NDMM, n = 56; RRMM, n = 115). Specific esterases exhibited relatively high or low expression in MM, and expression of specific esterases (UCHL5, SIAE, ESD, PAFAH1B3, PNPLA4 and PON1) was significantly altered on progression from NDMM to RRMM. High expression of OVCA2, PAFAH1B3, SIAE and USP4, and low expression of PCED1B, were identified as poor prognostic markers (P < 0.05). The MMRF CoMMpass dataset provided validation that higher expression of PAFAH1B3 and SIAE, and lower expression of PCED1B, were associated with poor prognosis. Conclusions Esterase gene expression levels change as patients progress from NDMM to RRMM. High expression of OVCA2, PAFAH1B3, USP4 and SIAE, and low expression of PCED1B, are poor prognostic markers in MM, suggesting a role for these esterases in myeloma biology.


BACKGROUND
Both antibody-drug conjugates (ADCs) and peptide-drug conjugates (PDCs) represent important therapeutic classes that enable the selective introduction of cytotoxic drugs into cancer cells over healthy cells, potentially improving efficacy and reducing systemic toxicity compared with non-conjugated versions of the same drug. 1,2 Conceptual similarities between ADCs and PDCs include selective targeting and subsequent cellular internalisation, although the exact mechanism of action is specific to each individual conjugate. To release cytotoxic payloads within the target cell, ADCs are typically endocytosed and cleaved within the consequent lysosomal structure, whereas PDCs undergo hydrolytic cleavage in the cell via interactions with cell surface receptors, or after direct entry into the cytosol through the cell membrane. [1][2][3] Various different enzyme classes may be involved in conjugate metabolism and release of the cytotoxic payload in ADCs and PDCs, and hence research interest in metabolising enzymes is growing for both therapeutic approaches. [1][2][3] The esterase enzyme family is a subclass of the hydrolase enzyme superfamily that functions to hydrolyse ester bonds. 4,5 Many different esterases have been identified, which differ in their substrate specificity and biological function, and the expression of specific esterases may be dysregulated in cancer. Hydrolysing enzymes can be highly expressed in cancer cells, and have previously been implicated in the reprogramming of metabolic pathways, promotion of cancer pathogenesis, drug metabolism and drug toxicity. 4,5 Esterases expressed in tumour cells may also differ in their stereoselectivity for hydrolysis of chiral esters compared with esterases expressed in healthy tissues. [6][7][8][9] Esterases can themselves also be administered to treat haematological malignancies; for example, the enzyme asparaginase has been used as an effective agent to treat acute lymphoblastic leukaemia for many years, 10 and it is also being investigated for the treatment of acute myeloid leukaemia. 11 Esterase hydrolysis, therefore, represents an interesting potential strategy for selective activation of anticancer drugs within cancer cells that overexpress esterases, or which express esterases with different specificities, while minimising toxic effects on healthy cells and tissues. 6 The study of esterase expression in cancer is still in its infancy, but interest is growing as a result of the development of novel PDCs such as melflufen (melphalan flufenamide), which utilises intracellular aminopeptidases and esterases to release a cytotoxic payload in multiple myeloma cells. 3,12 Several esterases have previously been reported to be dysregulated in specific cancers, with some also having the potential to be predictive or prognostic biomarkers. For example, low expression of neurexophilin and PC-esterase domain family member 4 (NXPE4) mRNA is a prognostic marker for shorter survival in patients with colorectal cancer, 13 and expression of platelet-activating factor acetylhydrolase 1B2 (PAFAH1B2) is inversely correlated with patient survival in pancreatic ductal adenocarcinoma. 14 In lung cancer, acetylcholinesterase (ACHE) may act as a tumour suppressor and is downregulated, 15 and high expression of granzyme A (GZMA) is significantly associated with an improved prognosis in multiple cancer types. 16 Sequencing studies have identified a number of genes that are mutated with high frequency in multiple myeloma, 17,18 although genes encoding esterases do not feature prominently among them. This is most likely because esterase genes do not have high-frequency mutation rates in multiple myeloma, which in turn suggests that regulation of esterase gene expression occurs through mechanisms other than mutations. Still, expression of esterase genes can be dysregulated in this disease: carboxylesterase 1 (CES1), CES2 and butyrylcholinesterase (BCHE) can both be overexpressed in multiple myeloma cells, 19,20 and BCHE is associated with a poor prognosis in multiple myeloma and dichotomously expressed, suggesting its tightly controlled regulation. 21 Also, expression of ubiquitin carboxyl-terminal esterase L1 (UCHL1) in samples from patients with newly diagnosed multiple myeloma (NDMM) correlated with a high-risk subgroup of myeloma. 22 Other studies have suggested possible mechanisms for the role of esterases in the pathogenesis of multiple myeloma; these include impaired oxidative/antioxidative balance due to reduced activity of serum paraoxonase-1 and arylesterase, 23 and inhibition of ubiquitinspecific protease 14 and ubiquitin-C-terminal hydrolase-5 (UCHL5). 24,25 In this study, we evaluated the biological significance of a panel of esterases in multiple myeloma using gene expression profiling, cytogenetic profiling and clinical outcome analyses. We demonstrate that specific esterases exhibit relatively high or low expression in multiple myeloma, that the esterase expression profile changes on the progression of the disease, and that high or low expression of individual esterases is associated with poor prognosis. Although single-nucleotide variants (SNVs) were rare in the esterase genes, several had duplication and deletion copy number alterations.

METHODS
Sample collection and plasma cell enrichment For gene expression profiling and evaluation of genomic variants in the Institute for Molecular Medicine Finland (FIMM) cohort, bone marrow aspirates were obtained from multiple myeloma patients after obtaining written informed consent and following protocols approved by an ethical committee of the Helsinki University Hospital Comprehensive Cancer Center, and in compliance with the Declaration of Helsinki. Matched patient skin biopsies were collected (also with informed consent) at the same time and from the same site as bone marrow aspirates, and in accordance with approved protocols, for constitutional DNA analysis. Bone marrow mononuclear cells were isolated by Ficoll-Paque gradient centrifugation (GE Healthcare), and CD138+ plasma cells enriched by immuno-magnetic bead selection (StemCell Technologies).
RNA sequencing and analysis RNA was extracted from CD138+ plasma cells using the AllPrep ® DNA/RNA/miRNA Universal or miRNeasy kits (Qiagen). RNA integrity was measured on an Agilent Bioanalyzer 2100 instrument; only samples with RNA integrity ≥7 were used for sequencing. Illumina-compatible RNA sequencing libraries were prepared using Scriptseq™ or Nextera technology and sequenced on Illumina HiSeq ® 1500 or 2500 instruments (Illumina). After preprocessing, filtered reads were aligned to the GRCh38 human reference genome using the STAR aligner tool. 26 Gene read counts were normalised using the Reads Per Kilobase of transcript per million mapped reads (RPKM) method. In total, 51 annotated esterase genes (Supplementary Table S1) were identified in the human genome (assembly GRCh38) utilising the Ensembl release 99 27 and NCBI 28 databases, using the search term 'esterase' and further confirming the molecular function (gene ontology) of identified genes. A cut-off value of >1 RPKM was used to filter expressed esterase genes. DEseq2 29 was used to identify variation in gene expression in newly diagnosed versus relapsed/refractory samples. The contribution of esterase gene expression to survival outcome was estimated by Kaplan-Meier analysis and performed using expression-based filtering. For each esterase gene, samples were stratified into 'high' (≥median expression) and 'low' (<median expression) expression groups. The significance of the difference between the two groups (high versus low expression) was deduced using a Mantel-Cox log-rank test.
Exome sequencing and cytogenetics The DNeasy ® Blood & Tissue kit or AllPrep ® DNA/RNA/miRNA Universal kit (Qiagen) was used to isolate genomic DNA from skin biopsies and CD138+ cells. The SeqCap ® EZ MedExome kit (Roche NimbleGen), SureSelect Clinical Research Exome kit or SureSelect Human All Exon V5 kit (Agilent Technologies) was used for exome capture. Sequencing was performed on HiSeq ® 1500 or 2500 instruments. VarScan2 somatic algorithm 30 was implemented for calling somatic mutations, and mutation annotations were performed using SnpEff 4.04 31 as described previously. 32 Gene copy number variants (CNVs) were identified using the CopyCat tool (https://github.com/chrisamiller/copycat). Cytogenetics data were generated using fluorescence in situ hybridisation technology as described previously, 33

Patient population
In total, 134 patients provided 171 bone marrow aspirates for the FIMM dataset: 56 samples for the NDMM subgroup and 115 samples for the relapsed/refractory multiple myeloma (RRMM) subgroup (see Supplementary Fig. S1 for a patient/sample flow chart). Individual patients could provide samples only at diagnosis (n = 49), at both diagnosis and relapse (n = 7), only at first relapse (n = 59) or at recurring relapses, i.e., multiple relapse samples from one patient (second relapse, n = 13; third relapse, n = 3; fourth relapse, n = 2; sixth relapse, n = 1). The median age was similar in the NDMM and RRMM subgroups, and there was a higher proportion of males (78/134, 58.2%) in the total population (Table 1). More patients in the RRMM subgroup (n = 20, 25.6%) had a 17p deletion compared with the NDMM subgroup (n = 5, 8.9%; Fisher exact test P = 0.015). Chromosome 1q gain was also more common in the RRMM subgroup (n = 44, 56.4%) compared with the NDMM group (n = 14, 25.0%; Fisher exact test P = 0.0004). In the RRMM subgroup, previous treatments included alkylating agents in 97.4%, bortezomib in 88.5% and immunomodulatory drugs in 79.5% of patients.
Esterase gene expression profile in multiple myeloma samples RNA extracted from 123 of 171 CD138+ plasma cell patient samples was suitable for RNA sequencing analysis, which included 41 samples from patients with NDMM, and 82 samples from patients with RRMM.
Esterase gene expression levels were ranked based on their abundance ( Fig. 1a and Supplementary Table S1). The most abundant esterase mRNAs were ovarian tumour suppressor range of −10.7 (ASPG) to −5.1 (NLGN1; Fig. 1a). Similar esterase mRNA expression patterns were observed (r = 0.9, P value = 2.2e-16) utilising the MMRF CoMMpass dataset (N = 892) (Supplementary Figs. S2, S3). In our FIMM dataset, the esterase genes clustered in four subgroups based on expression level, with group I having the highest and group IV having the lowest level of expression (Fig. 1b).
The mRNA expression profiles of several esterases were significantly different in samples from newly diagnosed patients versus those whose disease had relapsed or were refractory to treatment (Fig. 2). The expression levels of UCHL5, sialic acid acetyl esterase (SIAE), esterase D (ESD), PAFAH1B3 and PON1 were significantly higher (P ≤ 0.01; adjusted P < 0.1) in RRMM versus NDMM samples, and the expression level of PNPLA4 was significantly lower (P = 0.007). In the MMRF CoMMpass dataset using paired samples (NDMM, n = 39; RRMM, n = 45), we observed similar expression patterns, with the median expression of the majority of these esterase genes being higher in RRMM samples. Whereas gene PNPLA4 had a conflicting expression pattern with expression being higher in RRMM samples in the CoMMpass dataset. However, none of the genes were found to have P values ≤ 0.05 ( Supplementary Fig. S4) in the DEseq2 analysis. In addition, the expression of the esterase gene phospholipase A2 Group VII (PLA2G7; P < 0.001) was significantly different in NDMM versus RRMM samples using paired samples from the MMRF CoMMpass data (data not shown).
DEseq2 analysis also revealed that out of all the esterases (n = 6) predicted to be differentially regulated in NDMM versus RRMM samples in our dataset, only UCHL5 was upregulated in both RRMM samples as well as samples with 1q gain ( Supplementary  Fig. S5).
Prognostic significance of esterase expression in multiple myeloma In our FIMM dataset, patient samples exhibiting high expression of OVCA2, PAFAH1B3, SIAE and USP4 were associated with a significantly poorer prognosis versus those with low expression (Fig. 3a Supplementary Fig. S8).
Patient samples exhibiting low expression of GZMA, PCED1B and NXPE3 were associated with poorer prognosis versus those with high expression (Fig. 3a, Supplementary Figs Fig. S6).
In contrast, multiple esterase genes were found to have copy number alterations of both duplication and deletion types (Fig. 4). PNPLA6 was found to have a copy number gain in 30.8% (52/169) of samples, PAFAH1B3 had a copy number gain in 18.9% (32/169) of samples and CEL had a copy number gain in 17.2% (29/169) of samples. Deletions in genes encoding ESD, UCHL3 and ABHD13 were observed in more than 38% of the samples (Fig. 4a). These results were validated in the MMRF CoMMpass dataset: PNPLA6 was found to have a copy number gain in 40.1% (419/1044) of samples, CEL was found to have a copy number gain in 38.5% (402/1044) of samples and PAFAH1B3 had a copy number gain in 31.7% (331/1044) of samples. Deletions in genes encoding ESD, UCHL3 and ABHD13 were observed in~36-41% of samples ( Supplementary Fig. S11). Furthermore, in the FIMM cohort, a   correlation analysis of gene expression and CNVs revealed that among all genes predicted to have CNVs, only the genes UCHL5 and UCHL3 had a weak positive correlation (r > 0.5) of 0.57 and 0.54 with corresponding log2(RPKM) gene expression values (Fig. 4b).

DISCUSSION
A detailed understanding of the role of specific metabolic enzymes in the processing of novel therapeutic approaches such as ADCs and PDCs may be critical for their success in the treatment of specific cancers. The majority of chemical linkers in ADCs are based on hydrazone, disulfide, thioester or peptide bonds; these linkers are designed to exploit differences in intracellular pH, intracellular reduction potential or intracellular metabolic enzyme concentrations to break the linker and release a cytotoxic payload into tumour cells. 1 The cytotoxic activity of PDCs is also dependent on intracellular metabolic enzyme concentrations. 2 Therefore, among other enzymes, esterases have a significant potential for exploitation in drug design. However, their role in cancer, and more specifically in multiple myeloma, has not been widely studied to date, despite myeloma being a disease of altered protein homoeostasis that is commonly targeted with proteasome inhibitors. Novel lipophilic PDCs such as melflufen (melphalan flufenamide) can readily diffuse into multiple myeloma cells, where high expression of aminopeptidases results in hydrolytic cleavage and rapid release of the cytotoxic alkylator payload, and esterases can also create active intermediate metabolites. 3 Therefore, the role of these enzymes in myeloma cells is of particular interest; for example, aminopeptidase and esterase gene expression profiles could, hypothetically, be used to identify patient subgroups who are more likely to respond to drugs such as melflufen, which utilise these enzymes as part of their mechanism of action.
In this study, we demonstrate for the first time that several individual esterase genes exhibit relatively high or low median expression in bone marrow aspirates from patients with multiple myeloma. Interestingly, the expression profile of several genes appeared to change on progression from NDMM to RRMM; the expression of UCHL5, sialic acid acetyl esterase (SIAE), esterase D (ESD), PAFAH1B3 and PON1 was significantly higher in RRMM versus NDMM samples, whereas the expression of PNPLA4 was significantly lower in RRMM versus NDMM samples. Among these differentially expressed genes, only UCHL5 was shown to be upregulated in samples with 1q gain (a common cytogenetic abnormality in MM), which is likely due to the genomic location of UCHL5 (1q31.2). UCHL5 is a deubiquitylating enzyme (DUB) that is more highly expressed in MM cells than in normal plasma cells. 25 Preclinical studies have shown that inhibition of UCHL5 decreases viability and inhibits proliferation of MM cells, and overcomes resistance to bortezomib. 25 Many DUBs demonstrate esterase activity, and DUB inhibitors have been studied for anti-myeloma activity. One such inhibitor, VLX1570, has been assessed in a phase 1 study of patients with RRMM; 35 this study was discontinued due to severe pulmonary toxicity. Nevertheless, further efforts to identify DUB inhibitors with a wider therapeutic index may be warranted given promising preclinical anti-tumour effects and activity in MM resistant to proteasome inhibitors. 25,35 In our dataset, high expression of OVCA2, PAFAH1B3, SIAE and USP4 was associated with a significantly poorer prognosis compared with samples showing low expression, whereas low expression of GZMA, PCED1B and NXPE3 was associated with a significantly poorer prognosis compared with samples expressing higher levels of the enzyme. Esterase expression patterns were similar in the MMRF CoMMpass validation dataset, but of the genes identified as being prognostic in our dataset, only high expression of PAFAH1B3 and SIAE, and low expression of PCED1B was associated with poor prognosis in the CoMMpass dataset. Other esterases predicted to have a role in disease prognosis only in the validation dataset included NXPE4, UCHL5, PAFAH1B2, BPHL, NXPE1 and UCHL3. The differences observed between datasets may simply have been caused by the lower number of samples in our dataset, or by technical differences in RNA sequencing library preparation methods. Another possible reason for the differences may have been that our dataset includes samples from patients with NDMM (32.7%) or RRMM (67.3%), whereas the CoMMpass dataset includes mostly NDMM samples (87.4%). Despite these differences, we were still able to identify and validate three prognostic genes.
The other esterases identified as being prognostic in our dataset, but not in the validation dataset, have all previously been reported to be dysregulated in cancer. OVCA2 was originally identified as a tumour suppressor gene expressed in healthy surface epithelial cells of the ovary, and is downregulated in the majority of ovarian cell lines and tumours. 36 OVCA2 is located at chromosome 17p13, which is commonly deleted in haematological malignancies and also harbours several tumour suppressor genes, including TP53. In multiple myeloma, del17p is widely accepted to be one of the most aggressive features of the disease, and a marker of poor prognosis. [36][37][38][39] In our dataset, OVCA2 was more highly expressed in RRMM versus NDMM samples, whereas in the MMRF CoMMpass dataset there was no difference in expression level. One possible explanation for this difference could be that our dataset had a lower frequency of del17p in RRMM samples (19/82, 23.2%) compared to the MMRF CoMMpass dataset (14/45, 31.1%). Also, in contrast with del17p being a marker of poor prognosis, our analysis of the FIMM dataset suggested that high expression of OVCA2 is associated with poor prognosis in multiple myeloma, despite low expression of OVCA2 being associated with del17p. Further research is required to determine the role of OVCA2 in multiple myeloma. Upregulation of PAFAH1B3 expression in multiple cancers has been observed previously, 4,40,41 with selective inhibition impairing tumour cell survival. 42 Very little has been published regarding SIAE in cancer cells. A recent in vitro study investigating the potential of targeting the ganglioside GD3 acetylation pathway to treat medulloblastoma showed that SIAE may be associated with mitochondria-mediated apoptosis and etoposide sensitivity. 43 There are also very few publications on PCED1B. The role of PCED1B antisense RNA was investigated in an in vitro study of glioblastoma, which indicated that it functioned in a hypoxia-inducible factor (HIF)-1-dependent manner and has potential as a prognostic biomarker and druggable target for GBM. 44 Further investigation into the role of these esterases in multiple myeloma may be warranted. Our analysis did not take account of gene heterogeneity in CD138+ plasma cells, and the contribution of this heterogeneity to the observed differences in esterase gene expression. Further studies using single-cell RNA sequencing analysis would be useful in this regard. In RRMM patients, it is possible that prior treatment types also influenced the observed esterase gene expression patterns; this is again worthy of further study.
Interestingly, the frequency of SNVs appeared to be rare in the esterase genes, in both our dataset and the MMRF CoMMpass dataset. It has been previously reported that, except for a subset of specific genes, recurrent mutation rates are low in patients with multiple myeloma, suggesting that the dysregulation of key signalling pathways, rather than single-gene mutations, is the key driver for malignancy. 18,45 Also, multiple myeloma is thought to be genetically heterogeneous, whereby clonal diversity would result in specific mutations only being present in a small number of cells within a tumour, confounding the molecular characterisation of tissue samples. 18 However, in our study, several esterase genes were found to have a high frequency of duplication and deletion copy number alterations. To explore the possibility that tumours lacking SNVs might be driven by copy number alterations or chromosomal rearrangements, one study reported that of 153 patients with multiple myeloma, 119 (77.8%) had evidence of at least one focal gene copy number gain or loss within a significant peak, including 40 of 60 patients (66.6%) lacking somatic SNVs in the most significantly mutated genes. 18 Further studies are required to establish the significance of esterase gene copy number alterations in patients with multiple myeloma.
One limitation of our expression analysis was the sample size; 123 samples were collected from patients with multiple myeloma, preventing statistically robust comparisons with the RRMM or NDMM samples. Another limitation was that matched samples were not widely available for patients progressing from NDMM to RRMM. Finally, regarding verification of the significance of esterase gene expression using the MMRF CoMMpass registry, it should be noted that this dataset was generated from diagnostic samples rather than treatment samples, and that the majority of the samples were taken from patients with NDMM; the possibility that some of our findings are more relevant to RRMM cannot be excluded.
In conclusion, specific esterases exhibited relatively high or low expression in multiple myeloma, and the esterase gene expression profile appeared to change on progression from NDMM to RRMM. High expression of PAFAH1B3 and SIAE, and low expression of PCED1B were identified as poor prognostic markers in patients with multiple myeloma, suggesting a role for these esterases in myeloma biology. Further work is needed to elucidate the biological significance of esterases in cancer, to better understand how they can be effectively utilised to activate anticancer drugs in tumour cells, and their potential applicability as biomarkers of MM and its progression.

ADDITIONAL INFORMATION
Ethics approval and consent to participate The procedures in this study were approved by an ethical committee of the Helsinki University Hospital Comprehensive Cancer Center (study permit 303/13/03/01/2011) and performed in accordance with the Helsinki Declaration of 1975, as revised in 2008. All patients provided their written informed consent.
Consent to publish Not applicable.
Data availability All data generated or analysed during this study are included in this article (and its supplementary information files).  Correlation comparison of gene expression and CNV scores for the genes UCHL3 and UCHL5. CNV copy number variation. a A CNV score of more than 0.5 predicts a duplication/gain event and a CNV score of less than −0.6 predicts a deletion event. manuscript development. J.L. has received personal fees from Amgen, Bristol-Myers Squibb, Celgene, Janssen, Novartis, Pfizer, Roche and Takeda. R.S. has received research funding and personal fees from Amgen, Bristol-Myers Squibb, Celgene and Takeda, and personal fees from Sanofi. N.N.N. is a consultant for Oncopeptides AB. F.L. is an employee of Oncopeptides AB. C.A.H. has received research funding from Celgene, Novartis, Oncopeptides AB, Orion Pharma, Pfizer and the IMI2 consortium project HARMONY. P.A. has received personal fees from Amgen, Bayer, Bristol-Meyers-Squibb, Celgene, GSK, Janssen, Novartis, Sanofi and Takeda. R.K. and M.M.M. have nothing to disclose.