Gene expression profile identifies distinct molecular subtypes and potential therapeutic genes in Merkel cell carcinoma.

Merkel cell carcinoma (MCC) is a rare primary cutaneous neoplasm of neuroendocrine carcinoma of the skin. About 80% of the MCC occurs due to Merkel cell polyomavirus (MCPyV) and 20% of the tumors usually occur due to severe UV exposure which is a more aggressive type of MCC. It tends to have an increased incidence rate among elderly and immunosuppressed individuals. On therapeutic level, sub-classification of MCC through molecular subtyping has emerged as a promising technique for MCC prognosis. In current study, two consistent distinct molecular subtypes of MCCs were identified using gene expression profiling data. Subtypes I MCCs were associated with spliceosome, DNA replication and cellular pathways. On the other hand, genes overexpressed in subtype II were found active in TNF signalling pathway and MAPK signalling pathway. We proposed different therapeutic targets based on subtype specificity, such as PTCH1, CDKN2A, AURKA in case of subtype I and MCL1, FGFR2 for subtype II. Such findings may provide fruitful knowledge to understand the intrinsic subtypes of MCCs and the pathways involved in distinct subtype oncogenesis, and will further advance the knowledge in developing a specific therapeutic strategy for these MCC subtypes.

Introduction

Merkel cell carcinoma (MCC), though with a low occurrence rate is a lethal form of skin cancer. Merkel cell carcinoma resembles Merkel cells, somatosensory cells present in epidermis's basal layer. Merkel cells are believed to act as slow mechanoreceptors [1]. MCC has become the centre of attention due to the discovery of Merkel cell polyomavirus and rapid increase in occurrence rate [tripled since 1986] [2]. MCC's incident rate is higher in older people and occurs to those parts of the skin which are exposed to the sun, suggesting its occurrence related to UV exposure [3,4]. Oxygen radicals produced due to UV-A exposure cause DNA damage by chromosomal breakage resulting in loss or gain of chromosomal material [5]. MCC, though fully understood yet, is also thought to contribute towards the development of other malignancies such as malignant melanoma and chronic lymphocytic leukemia etc. [6,7] Patients with MCC developed in neck and head, and those with positive lymph nodes are initially subjected to local excision, biopsy of sentinel nodes and are later provided radiotherapy [8]. But in case of MCC metastasis, chemotherapy such as cytotoxic chemotherapy, is found to have less efficacy in most of the cases. However, the major drawback of chemotherapeutics is the high reoccurrence rate of the disease. Thus, various alternative therapies are searched to overcome the shortfalls of these traditional therapeutics. One of the alternatives to the traditional therapies is immuno-therapeutics, in which monoclonal antibodies are used to treat stage IV MCC [9,10]. Therefore, molecular target therapies are considered a viable option for the immunocompromised patients with failed track record against immunotherapies. Molecular subtyping based on gene expression patterns has allowed researchers to develop tools for the identification of molecular events associated subtypes in various tumors such as colon, lung, ovarian and prostate cancer, and provided somewhat success in prognosis and treatment [[11], [12], [13], [14], [15], [16], [17], [18]]. To date, limited data is available to correlate MCC with distinct molecular subtypes. In current study, two distinct molecular subtypes were identified through gene expression profile data. Various therapeutic genes and pathways were identified in the sample molecular subtypes which can be helpful in unfolding the novel target therapies for MCCs based on MCC molecular subtypes.

Material and methods

Determination of subtypes in MCC

Gene expression Omnibus (GEO) was used to collect the Expression profiling data of Merkel cell carcinoma cases. Dataset of GSE39612 which contains 30 MCC cases [ 19] and the second Dataset of GSE22396 having 35 MCC cases [20] were used to determine the molecular subtypes of MCC. For the identification of MCC molecular subtypes, each dataset was filtered with standard deviation and subsequently transformed on the basis of gene-based centring, then transformed dataset was separately ran on R package of Consensus Clustering Plus [21], having distance factors of 1-pearson correlation and 80% gene resampling, maximum evaluated k of 12, and agglomerative hierarchical clustering algorithm over 1000 iterations. The accuracy of Consensus Clustering Plus was checked by calculation of silhouette width (R package clustering) [22].

Reproducibility measurement of the MCC subtypes

Reproducibility of MCC molecular subtypes between the cohorts of GSE39612 and GSE22396, was determined through Subclass Mapping (SubMap) of the Gene Pattern. Subclass Mapping was carried out with parameters of (num. marker. genes = 300, num. perm = 1000 and num. perm. fisher = 1000) [23].

Gene ontology and gene set enrichment (GSEA) analysis

For the examination of Subtype specific genes, SAM [24] was performed with false discovery rate (FDR) of 0.05. The David Bioinformatics Resources online version 6.7 (https://david.ncifcrf.gov/) was preferred to emulate the GO and KEGG pathways. Whereas GSEA [25] was used to examined gene expression patterns and pathways of each subtypes. Potential therapeutic genes of each subtype of MCCs were explored through Target V2 database (http://www.broadinsitute.org/ccancer/cga/target).

Statistical analysis

Kaplan-Meier plots were plotted to calculate the survival curve. Whereas Chi-square and Fisher's exact tests were performed to evaluate the statistical significance of association between subtypes and clinical factors, and those having p value lower than 0.05 was considered as significant.

Results

Consensus clustering identifies two distinct MCC molecular subtypes

MCC molecular subtypes were determined through Consensus clustering by checking the genes with high variety of expression levels in MCC cohort. Two molecular subtypes were identified in GSE39612 by choosing the optimal number of clusters with highest increase in the area under curve of empirical cumulative distribution (CDF) (Fig. 1, A-C). Silhouette width determination was used to confirm the subtype assignment confidence and those found with positive silhouette value were selected and subsequently analyzed. Out of 30 MCC cases, 29 cases (core samples) were found to have positive silhouette values, 21 out of these 29 cases belong to subtype I and subtype II related cases were found to be 9 in number (Fig. 1D).

Fig. 1

Identification of two molecular subtypes of MCC in GSE39612. (A) Empirical cumulative distribution plot determines the optimal number of MCC molecular subtypes. (B) Relative increase in the area under the CDF curve along with increasing assumed number of molecular subtypes. (C) Consensus clustering matrix of MCC samples using two molecular subtypes. (D) Silhouette analysis of MCC samples based on the assignment from Consensus Clustering.

Validation of subtypes in the independent dataset

In order to confirm the validity of MCC molecular subtypes for clinic use, additional independent dataset (GSE22396) containing 35 cases was obtained from GEO as an independent validation cohort. As a result, two different molecular subtypes were observed through Consensus Clustering analysis as well (Fig. 2). Silhouette analysis revealed that 31 out of 35 samples had positive silhouette values and were defined as core samples. In 31 silhouette positive samples 21 belong to subtype I and 10 belong to subtype II.

Fig. 2

Identification of two molecular subtypes of MCC in GSE22396. (A) Empirical cumulative distribution plot determines the optimal number of MCC molecular subtypes. (B) Relative increase in the area under the CDF curve along with increasing assumed number of molecular subtypes. (C) Consensus clustering matrix of MCC samples using two molecular subtypes. (D) Silhouette analysis of MCC samples based on the assignment from Consensus Clustering.

Reproducible molecular subtypes in the two cohorts

SubMap analysis was used to measure whether molecular subtypes from the two datasets can be recapitulated or not. Results suggested the significant correlation between A1-A2 subtypes and B1-B2 subtypes of GSE22396 and GSE39612 respectively, indicating the presence of two consistent different MCC molecular subtypes with distinct gene expression patterns (Fig. 3).

Fig. 3

Association in the SubMap matrix between the subtypes of two independent dataset GSE3916 and GSE22396 showing the significant correlation. The correlation significance was denoted by FDR-corrected p-value.

Molecular subtypes and their relation with clinical symptoms

The relationship between molecular subtypes and clinical characteristics was analyzed based on GSE39612 (Supplementary Table S1). MCC molecular subtype I presented significantly high rate in upper site (Head and Neck) tumor (10/30, 33.33%) as compared to MCC subtype II (1/30, 3%) (p = 0.0468⁎). In this dataset of MCC, the number of female patients (46.6%) is higher than the male (43.3%) patients. Between the two molecular subtypes a significant difference was observed among the patients on the basis of their age group (p = 0.0273⁎). The mean age of subtype I was 79 years old which was higher than subtype II patients which was 70 years old. The median time of overall survival (OS) for subtype I MCCs was found to be 13.5 months which was lower than the patients with subtype II MCCs who had 16.5 months, although the difference is not significant (p = 0.507). In MCC subtype I, 38% of the patients were found to be MCPyV positive while the 57.1% were MCPyV negative. Whereas in subtype II of MCC 25% of the patients were MCPyV positive and 50% were MCPyV negative (p = 0.483). In addition, no remarkable difference was observed in the Kaplan-Meier plots (survival curve analysis) (Supplementary Table S1).

Functional enrichment analysis of MCC subtype-specific genes

By SAM analysis of GSE39612, 28,912 genes were found to be significantly differentially expressed between the two MCC molecular subtypes. According to the results, 10,920 out of 28,912 genes were specifically overexpressed in subtype I MCCs as compared to subtype II MCCs, whereas, the remaining 17,992 genes were found overexpressed in subtype II MCCs (Supplementary Table S2). GO and KEGG analysis were performed to obtain the biological information of the subtypes. Top 1000 most overexpressed genes in each MCC subtype were selected. Out of total, 145 biological processes and 19 KEGG pathways were found to be enriched in subtype I MCCs, whereas 139 biological processes and 31 KEGG pathways were found to be enriched in subtype II MCCs. Remarkably, 43 (4.3%) genes of the subtype I MCC specific top 1000 genes were found to be related to DNA replication pathway (Supplementary Table S3). Additional pathways found enriched in subtype I MCCs were spliceosome, Mismatch repair and Cell cycle etc. (Fig. 4A). On the other hand, genes overexpressed in subtype II MCCs were related with distinct biological pathways such as TNF signalling pathway, MAPK signalling, myeloma formation and other regulatory pathways etc. (Fig. 4B).

Fig. 4

Pathways enriched in each MCC subtypes. (A) KEGG pathways in subtype I. (B) KEGG pathways in subtype II.

Potential clinical implication of MCC subtyping

The basic aim behind the determination of MCC molecular subtypes is to develop various therapeutic routes specific towards distinct subtypes and their further use in clinical studies and discourses. For the determination of therapeutic modules related to each MCC subtype, overexpressed genes based on subtype specificity were compared with target Database (http://www.braodinsitute.org/cancer/cga/target) containing gene targets and functional inhibitors [26], for further studies to be carried out on the gene targets to convert them into potential clinical targets [[27], [28], [29]].

In total, 11 genes were found to be specific against molecular subtypes and would benefit subtype-specific MCC patients (Table 1). Subtype I MCCs contains PTCH1, SMARCA4, CDKN2A, AURKA and BRCA1, whereas subtype II MCCs contains 6 overexpressed genes namely MCL1, FGFR2, BRD4, MET, SMO and PDGFRA.

Table 1

Gene overexpressed in each MCC molecular subtype.

Gene overexpressed		Examples of potential therapeutic agents
Subtype I	PTCH1	Vismodegib, hedgehog inhibitors
	SMARCA4	HDAC
	CDKN2A	CDK4/6 inhibitors
	AURKA	AURKA inhibitors
	BRCA1	PARP inhibitor
	MCL1	Tubulins
	FGFR2	FGFR inhibitors
Subtype II	BRD4	HDAC inhibitorsBromodomain inhibitors
	MET	Gefitinib, erlotinib, EGFR inhibitorsCrizotinib, MET inhibitors
	SMO	Vismodegib, hedgehog inhibitors
	PDGFRA	Imatinib

Discussion

MCC is an aggressive type of tumor found in the skin with higher metastasis, recurrence, and mortality rate. Since last two decades, the incidences of MCC have nearly tripled. Overall, 5-year survival rate of MCC patients is almost 30–60% [30]. Advanced age, UV exposure and immunosuppression increase the risk of MCC [31]. Usual course of treatment is radiation therapy followed by wide surgical excision for most of the cases, but in-case of metastasis and positive lymph nodes, only chemotherapy is used [[32], [33], [34]].

The survival rate of MCC patients is extended if it is diagnosed early. Molecular subtyping of tumors based on gene expression profiling has guided subtype-specific diagnosis, prognosis, and aided in developing subtype targeted therapies [21]. The clinical trial of Herceptin treatment in breast cancer is an ample example of subtype targeted therapies. Wherein, HER2-negative breast cancer patients did not get any benefit from Herceptin treatment yet HER2-positive breast cancer patients responded very well and benefited from this treatment [35]. Therefore, analysing the molecular heterogeneity of MCC is posed to provide deep knowledge of genes and pathways involved in MCC. It will also provide the new opportunities to target subtype specific patients of MCC.

In our study, we have identified two distinct MCC molecular subtypes. Gene ontology and gene set enrichment analysis revealed distinct gene signatures in both subtypes of MCC. Spliceosome, Cell cycle, DNA replication and Mismatch repair etc., were the pathways over-expressed in subtype I MCCs (Fig. 4) along with the overexpression of DEPDC1, ELAVL and OIP5 genes in subtype I (Fig. 5). DEPDC1 is a cell cycle related gene in the bladder cancer [[36], [37], [38], [39]]. Dysregulation of DEPDC1 plays a significant role in the regulation of cell cycle, motility and progression of nasopharyngeal carcinoma [[40], [41], [42]], as well as in other human cancers (colorectal cancer, hepatocellular carcinoma and glioblastoma) [[43], [44], [45]]. Whereas, ELAVL1 gene encodes HuR protein. The important function of HuR is for mRNA stability. High levels of cytoplasmic HuR were observed in numbers of various cancers [[46], [47], [48], [49]]. For instance, in glioma, it was found involved in tumor growth and in the onset of drug resistance [[50], [51], [52], [53], [54], [55], [56], [57]]. A protein coding gene, Opa interacting protein 5 (OIP5) belongs to cancer/testis antigens (CTAs) [58]. OIP5 is found involved in different biological processes of tumor [59], Gong et al. found overexpression of OIP5 gene in clear renal cell carcinoma [60]. On the basis of above gene expression and cellular pathways found in subtype I MCCs of this cohort, potential therapies based on established roles of such determinants can be developed for subtype I MCC patients. In contrast, gene ontology and gene set enrichment analysis revealed the overexpressed pathways found in subtype II include TNF signalling pathway, MAPK signalling pathway, Pathways in cancer (Supplementary Table S4). The expression, repression and re-expression cycle of genes like H19 (Human Long Noncoding) was found to be an important oncogenic factor in several cancer tissues such as bladder cancer and hepatocellular carcinoma etc. [[61], [62], [63], [64], [65]] Whereas, expression pattern of RGS1 gene plays an important role in pathogenesis of the different malignancies and regulation of chemokine induced signalling of B cells [[66], [67], [68], [69], [70], [71]]. Overexpression of these genes in subtype II MCCs and their roles in other cancers may open new doors of opportunities to understand the patients of subtype II MCC.

Fig. 5

GSEA reveals different gene expression signature in distinct MCC molecular subtypes. (A) Representing different gene expression patterns in subtype I and subtype II. Red, overexpressed genes; blue, down expressed genes. (B) GSEA shows the activity of DNA replication and Spliceosome pathways in subtype I. (C) GSEA demonstrated the activity of Toll like receptor signalling and cytosolic DNA sensing pathways in subtype II. Whereas NES is denoting normalized enriched score and FDR is denoting false discovery rate.

We also identified a number of known target genes and their potential therapeutic agents using TARGET V2 database in subtype I and II MCCs. These genes include PTCH1, CDKN2A, AURKA and BRCA1 in subtype I, and MCL1 in subtype II (Table 1). Cellular pathways such as Hh pathway is involved in the formation of organs, tissues and motor neurons [72]. Under biological conditions, the Hh signalling pathway is regulated through a transmembrane receptor called Patched 1 (PTCH1). Hh receptor Ptch1 is found to be overexpressed in many types of cancers such as Breast, prostate, melanoma, lung cancer and myeloid leukemia [[73], [74], [75], [76], [77]]. Interestingly, as a novel small molecule inhibitor of Hh, vismodegib was found actively involved in generation of anti-tumor response in advanced Basal cell carcinoma and got approval by the US FDA priority review program on January 30, 2012 [[78], [79], [80], [81]]. However, in another study, Carroll et al. suggested that Hh inhibition could not be an active therapy for MCC [82]. In current study, we found that more than 60% patients of both datasets GSE39612 and GSE22396 belong to subtype I of MCC (Fig. 1D, Fig. 2D and Supplementary Table S1). Our results revealed the importance of PTCH1 inhibitors in MCC subtype I (Table 1). Similar outcome was also reported in two other studies indicating the importance of Hh inhibitors in MCC. They observed the presence of amplified Hh pathway components in MCC tumors and proposed the importance of Hh pathway in pathogenesis of MCC [83]. Which further warrant the use of Hh inhibitors in the treatment of MCC [84]. We have also identified Histone deacetylase inhibitors (HDACi) in patients of subtype I MCC. Function of HDACis is to upregulate the low HLA class-I expression in MCC patients to contain the resistance of MCC towards PD-1/PD-L1 inhibition and low influx of CD8+ T which allows the tumor to avoid host immune response [85]. Thus, possible use of PTCH1 and HDACi inhibitors might be beneficial for the patients of subtype I MCC. Moreover, extension of antiapoptotic BCL2 family proteins (MCL1, BCL2, and BCL-xL) has been reported in relation to chemotherapy resistance in various types of cancers [86]. MCL1 is an influencing therapeutic target in cancer. Various gene knockdown and pharmacological inhibition models have demonstrated the potential of MCL1 as a critical pro-survival factor in tumors such as AML, multiple myeloma, B-cell lymphomas, and breast cancer [[87], [88], [89], [90], [91], [92]]. Microtubules also play a fundamental role in various cellular processes and are favourite antitubulin target chemotherapies in different cancers. Microtubule inhibitory target agents such as Vinscristine and Taxol are widely used in antitubulin chemotherapeutics. They are able to cause mitotic arrest in cells and ultimately lead towards apoptosis. MCL1 being an important regulator in the process of apoptosis is a key target of antitubulin chemotherapies [93]. In one study, 10 out of 11 MCC xenografts, ABT-263 successfully inhibited the multiple pro-survival Bcl-2 proteins and strongly support the therapeutic effect of these molecules MCC targeting [94]. Another pan-active BCL2 inhibitor, Sabutoclax was found capable of targeting all 5 anti-apoptotic proteins. A1210477 and S63845 also followed up the promise that MCL1 inhibition seems to offer as MCL1 specific inhibitors [95]. According to our gene expression analysis, we found that MCL1 inhibitor played a significant role in subtype II MCC under the influence of BCL-2 family proteins. Similarly, one more prospective inhibitor against targeted MCC called Imatinib was also identified. Imatinib inhibitor therapy is considered as an efficient therapy for patients with MCC metastasis and complications of inoperability [96]. These results suggested that patients of Subtype II MCC might benefit from these inhibitors.

In conclusion, we defined distinct molecular subtypes of MCCs. Our findings provided subtype-specific mechanisms underlying tumorigenesis and tumor progression. These outcomes might help in better targeting of the disease and developing an individualized medicine through subtype-specific manner for MCCs patients.

Author contributions statement

Xiangqian Guo: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Resources, Data Curation, Writing - Original Draft, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. Umair Ali Khan Saddozai: Validation, Formal analysis, Investigation, Data Curation, Writing - Original Draft, Writing - Review & Editing, Visualization. Fengling Wang: Formal analysis. Zhang Lu: Formal analysis. Xinying Ji: Resources, Writing - Review & Editing, Supervision, Project administration. Yongqiang Li: Resources, Supervision, Project administration. Muhammad Usman Akbar: Writing - Original Draft, Writing - Review & Editing. Wan Zhu: Writing - Original Draft, Writing - Review & Editing. Yu Cheng: Writing - Review & Editing.

The following are the supplementary data related to this article.

Supplementary Table S1

Clinopathological characteristics of GSE39612 (N = 30).

Supplementary Table S2

Results of SAM analysis between different subtypes of MCC in GSE39612. Data is representing over expression of gene in two subtypes in positive and minus value.

(1). Gene were ordered according the SAM significance, gene with the positive value were overexpressed in subtype II of MCC whereas minus value genes were overexpressed in subtype I.

(2). Significance increased with decreased of absolute values of ranks from 10 to 1 or −10 to 1. Rank 1 or −1 genes are the most significant of SAM results.

Supplementary Table S3

Biological process enriched in each subtypes of MCC.

Supplementary Table S4

Enriched KEGG pathways in each subtypes of MCC.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.