Activation of Oncogenic and Immune-Response Pathways Is Linked to Disease-Specific Survival in Merkel Cell Carcinoma

Simple Summary Merkel cell carcinoma (MCC) is a rare and aggressive skin cancer. Developing targeted therapies for MCC requires increased understanding of the mechanisms driving tumor progression. In this study, we aimed to identify genes, signaling pathways, and processes that play crucial roles in determining disease-specific survival in MCC. We analyzed the gene expression of 102 MCC tumors and identified genes that were upregulated among survivors and in patients who died from MCC. We cross-referenced these genes with online databases to identify the pathways and processes in which they function. Genes upregulated among survivors were mostly immune response related and genes upregulated among patients who died from MCC function in various pathways that promote cancer progression. These results could guide future studies investigating whether these genes and pathways could be used as prognostic markers, as markers to guide therapy selection, or as targets of precision therapy in MCC. Abstract Background: Merkel cell carcinoma (MCC) is a rare but highly aggressive neuroendocrine carcinoma of the skin with a poor prognosis. Improving the prognosis of MCC by means of targeted therapies requires further understanding of the mechanisms that drive tumor progression. In this study, we aimed to identify the genes, processes, and pathways that play the most crucial roles in determining MCC outcomes. Methods: We investigated transcriptomes generated by RNA sequencing of formalin-fixed paraffin-embedded tissue samples of 102 MCC patients and identified the genes that were upregulated among survivors and in patients who died from MCC. We subsequently cross-referenced these genes with online databases to investigate the functions and pathways they represent. We further investigated differential gene expression based on viral status in patients who died from MCC. Results: We found several novel genes associated with MCC-specific survival. Genes upregulated in patients who died from MCC were most notably associated with angiogenesis and the PI3K-Akt and MAPK pathways; their expression predominantly had no association with viral status in patients who died from MCC. Genes upregulated among survivors were largely associated with antigen presentation and immune response. Conclusion: This outcome-based discrepancy in gene expression suggests that these pathways and processes likely play crucial roles in determining MCC outcomes.


Introduction
Merkel cell carcinoma (MCC) is a neuroendocrine carcinoma of the skin with a poor prognosis. The survival rate of MCC varies significantly; population-based studies from New Zealand, Finland, and the United States revealed 5-year disease-specific survival rates

Patients, Clinical Data, and Tissue Samples
Data on patients diagnosed with MCC in Finland from 1983 to 2013 were obtained from the Finnish Cancer Registry and Helsinki University Hospital files. Clinical details were extracted from hospital records. Formalin-fixed, paraffin-embedded (FFPE) tissue blocks from primary tumors were retrieved from the pathology archives. MCC diagnoses were confirmed in a blinded fashion from our earlier studies according to well-established criteria by two researchers with special expertise in MCC pathology [14]. The study protocol was approved by the Ethics Committee of Helsinki University Central Hospital. The Ministry of Health and Social Affairs granted permission to collect patient data and the National Authority for Medicolegal Affairs to collect tissue samples.
The patients were subdivided into the following groups: the poor prognosis group (patients who died from MCC) and the good prognosis group (patients who were either alive at the most recent follow-up or died from a cause unrelated to MCC). The genes that were upregulated in the tumors of patients in the poor prognosis group when compared with the good prognosis group will be referred to as death-associated genes (DAGs); genes that were upregulated in tumors of patients in the good prognosis group when compared with the poor prognosis group will be referred to as survival-associated genes (SAGs).
MCPyV detection from paraffinized tumor blocks was performed previously and is described in detail elsewhere [4]. Briefly, the presence of MCPyV DNA was analyzed from DNA extracted from representative deparaffinized tumor sections. Quantitation of MCPyV DNA was performed using real-time polymerase chain reaction (PCR). The relative DNA sequence copy number for each tissue sample was expressed as a ratio of MCPyV DNA to protein tyrosine phosphatase gamma receptor gene DNA. The sample was considered positive when MCPyV DNA copy number per reference gene was greater than 0.1.

RNA Extraction from FFPE Samples
RNA extraction and sequencing were performed on FFPE primary tumor samples from 120 patients from whom sufficient representative tumor tissue was available. Fourteen patients were excluded from the study as the corresponding samples did not pass the quality control during the processing of sequencing data. A further four patients were excluded because their tumor MCPyV status was not known. The final sample size was 102 patients.
Two 10 µM sections from FFPE MCC samples of a cancer-representative area were used as starting material for the RNA extractions. RNA extraction was performed with a QIASymphony SP instrument (QIAGEN GmbH, Hilden, Germany) and QIAsymphony RNA extraction kit (cat. No: 931636, QIAGEN GmbH, Hilden, Germany), following the manufacturer's protocol. The quality and quantity of the extracted RNA were measured using a 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA, USA). The average RNA integrity number (RIN) was 2.1.

3 RNA Sequencing
RNA sequencing was performed by the sequencing unit of the Institute for Molecular Medicine Finland (FIMM) Technology Centre, University of Helsinki. Prior to library preparation, 1 µL of 1:1000 ERCC RNA spike-in control was added to each sample. A QuantSeq 3 mRNA-Seq Library kit (Lexogen, Vienna, Austria, version 015UG009V0221) was used to prepare the RNA sequencing library in 48 sample batches. The library preparation was performed according to the manufacturer's instructions and is described in detail elsewhere [15]. The sample libraries were homogenous in average size and concentration. The sample libraries were pooled together by their molarity for the sequencing with 39 samples per lane. A HiSeq 2500 instrument from Illumina was used as the sequencer with high-output mode and v4 chemistry. The sequencing run was paired-end with a read length of 101 bp and one mismatch allowed in demultiplexing. The RNA sequence data were submitted to the SRA database in NCBI and can be accessed under the BioProject PRJNA775071.

Processing of Sequencing Data
Htseq-count files from Bluebee sequencing process were read into R [16] (version 4.1.0) and were matched to clinical data per sample. Gene annotation was retrieved from AnnotationHub (snapshotDate: 20 October 2021). Data were imported into DGEList object using edgeR package [17] (version 3.36.0). Genes were further filtered by filterByExpr function (default parameters) from edgeR. log2 counts per million were calculated with cpm function (default parameters) from edgeR. Additional sample quality control was performed by excluding samples with median logcount with deviation more than 10% from the median logcount over all samples. Samples without MCPyV status in clinical data were excluded as the last step of sample-based filtering.
In order to validate that our transcriptomic data derived from FFPE samples are comparable to data from fresh frozen tissues, we performed a comparative analysis between the data used in this study and data from Harms et al. (GSE39612) [18]. Data from Harms et al. were retrieved from GEO as GSE matrix with log2 RMA normalized signal intensity per sample per Affymetrix probeset. Intensities were median merged per gene symbol as per GPL data provided by GEO for the study in question. Only data from tumors with a known MCPyV status were included. We then performed a differential expression test of genes between MCPyV-positive and MCPyV-negative tumors separately in GSE39612 data and FFPE derived data. In the case of FFPE data, this was done with DESeq2 R package, and for GSE39612, it was done with limma R package. Of the top 250 (based on B-statistics) differentially expressed genes from the analysis of GSE39612, we further analyzed the 136 genes that were present in the FFPE dataset. LogFC-values of these 136 genes between MCPyV-positive and MCPyV-negative tumors in both datasets were plotted. The correlation between the logFC values of these 136 genes in the two datasets was calculated with the sm_statCorr() function in R with Pearson as the correlation metric.

Identification of DAGs and SAGs
A differential expression test of genes between the good and poor prognosis groups was performed using estimateDisp and exactTest function from edgeR with default parameters. DAGs were defined as those with a positive logFC with Benjamini-Hochberg (BH) corrected p-value < 0.05; these were genes found to be upregulated in samples of the poor prognosis group. SAGs had a negative logFC with BH corrected p-value < 0.05; these were genes found to be upregulated in samples of the good prognosis group. Altogether, we identified 79 genes differentially expressed between patients in these two groups.

Clustered Heatmap of Differentially Expressed Genes
To visualize the gene expression patterns of these 79 identified genes, we drew a conventional heatmap with 102 patients on the x-axis and the 79 genes on the y-axis. Both axes were clustered using 1-Pearson correlation coefficient as the distance measure and the linkage method Ward.D2. The rows of the matrix were scaled to have a mean of 0 and a standard deviation of 1.

Gene Ontology and Signaling Pathway Analyses of DAGs and SAGs
The gene ontology and signaling pathway analyses of DAGs and SAGs were performed by separately cross-referencing the lists of DAGs and SAGs with the Gene Ontology Biological Process (GO BP), the Gene Ontology Molecular Function (GO MF), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway databases using Enrichr [19,20]. Enrichr is a publicly available gene set search engine (http://amp.pharm.mssm.edu/Enrichr, accessed on 31 December 2021) [21]. The enriched GO terms and KEGG pathways were subsequently ranked by significance according to their p-values.

Statistical Analysis
Statistical analysis was performed with SPSS statistics 26.0 software (IBM Corporation, New York, NY, USA). The MCC-specific prognostic effects of tumor MCPyV status, sex, age at the time of diagnosis, tumor location, and the presence of lymph node or systemic metastasis at diagnosis were assessed by Kaplan-Meier analysis and logrank test. The factors that were significantly associated with MCC-specific survival in univariate analyses were included in a Cox regression multivariate analysis. In both uni-and multivariate survival analyses, MCC-specific survival was calculated from the date of diagnosis to death considered to be due to MCC, censoring subjects alive on their last follow-up date, and subjects who died from a cause unrelated to MCC. The correlation between the abundance of tumor-infiltrating, CD3-positive lymphocytes and the expression levels of the genes CD74, TAP2, PSME1, HLA-DRA, and HLA-E was investigated by calculating Pearson's product-moment correlation in each case. The abundance of tumor-infiltrating lymphocytes was investigated previously and is described in detail elsewhere; data were available in 80/102 cases [6].

Survival of ≥3 Years as Inclusion Criterion for the Good Prognosis Group
In order to evaluate the consistency of our findings, we repeated the differential expression test of genes between the good and poor prognosis group, this time only including the 40 patients who survived at least 3 years after initial diagnosis in the good prognosis group. We subsequently cross-referenced the differentially expressed genes yielded by this control analysis with the GO and KEGG databases in the manner described in Section 2.6.

Differential Expression Test of DAGs Based on MCPyV Status
In order to investigate potential relationships between DAG-expression and negative MCPyV status, we repeated the differential expression test of genes between MCPyVpositive (n = 15) and MCPyV-negative (n = 13) tumors of patients who died from MCC.

Overview of Patients
Detailed patient clinical data are shown in Table 1. MCC-specific death occurred in 28/102 (27%) cases and occurred within 5 years of initial MCC diagnosis in 25/28 (89%) cases. The median survival time before MCC-specific death was 2.3 years (range 0.3 to 14.8 years). There was one extreme outlier who survived 14.8 years before MCC-specific death; the second longest survival time was 6.3 years. The median follow-up time before death from a cause unrelated to MCC, or until the last follow-up date at which the patient was still alive, was 4.4 years (range 0.04 to 33.2 years). Kaplan-Meier analysis revealed that MCPyV-negativity, male sex, and the presence of lymph node or systemic metastasis at diagnosis correlated with poorer MCC-specific survival (all p < 0.001). Table 2 shows the results of the Cox regression multivariate analysis; all three variables were significantly associated with poorer MCC-specific survival in the multivariate analysis.

Identification of DAGs and SAGs
The differential expression test of genes revealed a total of 50 DAGs and 29 SAGs with a BH corrected p-value < 0.05. A summary of the upregulated genes, along with their corresponding logFC (fold change) values, false detection rates (FDR), and p-values, can be found in Table 3. A heatmap illustrating the expression of DAGs and SAGs across all samples can be found in Figure 1. Table 3. Genes associated with disease-specific death and survival in a series of 102 Merkel cell carcinoma patients.

Death-Associated Genes
Survival-Associated Genes

GO Enrichment and KEGG Signaling Pathway Analysis of DAGs and SAGs
The top 10 GO BP and GO MF terms most significantly enriched by the DAGs and SAGs, respectively, can be found arranged according to their p-values in Tables 4 and 5, along with their corresponding q-values and a list of the specific genes causing the enrichment. The q-value is an adjusted p-value calculated using the BH method for correction for multiple hypotheses testing. The top 10 KEGG pathway terms most significantly enriched by DAGs and SAGs can be found in the same format in Table 6.

GO Enrichment and KEGG Signaling Pathway Analysis of DAGs and SAGs
The top 10 GO BP and GO MF terms most significantly enriched by the DAGs and SAGs, respectively, can be found arranged according to their p-values in Tables 4 and 5, along with their corresponding q-values and a list of the specific genes causing the enrichment. The q-value is an adjusted p-value calculated using the BH method for correction for multiple hypotheses testing. The top 10 KEGG pathway terms most significantly enriched by DAGs and SAGs can be found in the same format in Table 6.
Cancer-relevant processes, functions, and pathways enriched by the DAGs included the PI3K-Akt signaling pathway, the MAPK signaling pathway, and angiogenic processes. Considering angiogenesis, the most significantly enriched GO terms included regulation of vascular associated smooth muscle cell migration, positive regulation of vascular endothelial cell proliferation, vascular endothelial growth factor receptor 2 binding, and vascular endothelial growth factor receptor binding. Among the most significantly enriched KEGG pathway terms was VEGF signaling pathway and among the DAGs was vascular endothelial growth factor A (VEGFA). Five of the DAGs were associated with the KEGG pathway term PI3K-Akt signaling pathway (MYB, AKT3, IGF2, ITGA6, and VEGFA) and four of the DAGs were associated with the KEGG pathway term MAPK signaling pathway (MECOM, AKT3, IGF2, and VEGFA).
The processes, functions, and pathways most significantly enriched by the SAGs were mostly related to immune response, particularly to antigen processing and MHC-dependent antigen presentation, and to T-cell mediated immune response. Examples include the KEGG pathway term antigen processing and presentation (CD74, TAP2, PSME1, HLA-DRA, HLA-E) and the GO BP terms T cell receptor signaling pathway (PSME1, HLA-DRA, LCP2, CARD11) and positive regulation of lymphocyte proliferation (SASH3, CD74, HLA-E).
Complete lists of all the GO BP, GO MF, and KEGG pathway terms that had p-values of <0.05, enriched by the DAGs and SAGs, are presented in Supplementary Tables S1 and S2. Table 4. Gene Ontology (GO) terms most significantly enriched by death-associated genes.  Table 5. Gene Ontology (GO) terms most significantly enriched by survival-associated genes.

Survival of ≥3 Years as Inclusion Criterion for the Good Prognosis Group
Including only patients who were still alive 3 years after diagnosis in the good prognosis group resulted in an increased amount of differentially expressed genes. The number of genes upregulated in the good prognosis group increased to 82, including 22/29 of the original SAGs; the number of genes upregulated in the poor prognosis group increased to 118, including 40/50 of the original DAGs. Considering the SAGs of this control analysis, the vast majority of the most significantly enriched GO and KEGG pathway terms were still immune response-related, especially related to antigen processing and MHC-dependent antigen presentation. The GO and KEGG pathway terms most significantly enriched by the DAGs in this control analysis still included embryonic developmental processes. The KEGG pathways PI3K-Akt and MAPK signaling pathway were still enriched by the same DAGs, with the exception of FGFR2 instead of VEGFA in both cases; however, the p-values were no longer significant owing to the increased number of differentially expressed genes. The DAGs in this control analysis also exhibited enrichment of several terms related to mitosis as well as chromatin binding, the latter in part caused by the inclusion of several genes encoding histone proteins; the main finding that was lost in this control analysis was the DAG VEGFA, and thereby several GO and KEGG pathway terms related to angiogenesis. Complete lists of the differentially expressed genes yielded by this control analysis, as well as the GO and KEGG pathway terms significantly enriched by them, are provided in Supplementary Tables S3-S5.

Differential Expression Test of DAGs Based on MCPyV Status
The differential expression test of genes between MCPyV-positive and MCPyV-negative tumors of patients who died from MCC yielded 100 genes significantly upregulated in MCPyV-negative tumors. Of these genes, eight (VSIG8, NEDD4L, ITGA6, KIF23, IGF2, H1-3, DST, COL21A1) were among the DAGs. A complete list of these 100 genes is provided in Supplementary Table S6.

Correlation between Tumor-Infiltrating Lymphocytes and SAGs Related to Antigen Processing and Presentation
The correlation analysis between an abundance of CD3-positive, tumor-infiltrating lymphocytes and the expression levels of the five SAGs causing the enrichment of the KEGG pathway antigen processing and presentation (CD74, TAP2, PSME1, HLA-DRA, and HLA-E) revealed a significant and positive Pearson's product-moment correlation in each case. Detailed results are presented in Supplementary Table S7.

Comparison of FFPE Data to Data from Fresh Frozen Tissues
Out of the 250 genes most significantly differentially expressed between MCPyVpositive and MCPyV-negative tumors in GSE39612, 136 genes were present in our FFPE data. The correlation analysis between the logFC-values of these 136 genes based of MCPyV status in GSE39612 and our data yielded a Pearson correlation coefficient of 0.82 (p < 0.001). A list of these 136 genes, together with their logFC-values in GSE39612 and our data, as well as a scatter plot of their expression, is provided in Supplementary Table S8.

Discussion
We found a dichotomous gene expression profile between the tumors of the poor and good survival groups. Many of the DAGs, which were upregulated in the poor survival group, function in various oncogenic pathways and processes. The SAGs, which were upregulated in the good survival group, were to a large extent immune-responserelated genes.
We found a clear association of DAGs with angiogenesis. Angiogenesis plays a crucial role in the progression of solid tumors, and increased tumor vascularization is a factor predicting poor prognosis in MCC [22]. VEGFA, a proangiogenic growth factor that was among the DAGs, has been found to be expressed in the majority of MCC tumors based on immunohistochemistry results. Overexpression of VEGFA has been shown to correlate with metastatic tumor spread of MCC [23,24].
Five of the DAGs were associated with the KEGG pathway term PI3K-Akt signaling pathway, and four of the DAGs were associated with the KEGG pathway term MAPK signaling pathway. Furthermore, although not recognized in the KEGG pathway database, the DAGs SULF2, RBFOX3, and COL11A1 have been reported to upregulate the PI3K-Akt signaling pathway and the DAGs ITGA6, SMYD3, and ENC1 have been reported to upregulate the MAPK signaling pathway in other malignancies [25][26][27][28][29][30][31]. MCPyV small tumor antigen has previously been demonstrated to activate p38 MAPK signaling in MCC, and immunohistochemical findings have demonstrated high degrees of activating AKT phosphorylation in MCC [32,33]. The PI3K-Akt and MAPK signaling pathways both serve oncogenic roles in several malignancies, such as stimulating cellular proliferation, inhibiting apoptosis, promoting tumor invasion and metastasis, and stimulating angiogenesis, which in the case of the MAPK pathway is partially mediated by upregulation of VEGFA [34,35].
Among the DAGs were also insulin-like growth factor 2 (IGF2), a protumorigenic growth factor, and IGFBP5, which can either inhibit or potentiate insulin-like growth factorsignaling depending on the context [36]. Insulin-like growth factor binding, insulin-like growth factor I binding, and insulin-like growth factor II binding were among the most enriched GO MF terms. MCC has previously been shown to express insulin-like growth factor-I receptor, but to the best of our knowledge, insulin-like growth factor signaling in MCC has not otherwise been reported [37].
Considering clinical implications, the aforementioned pathways and processes include several potentially viable targets for pharmacological intervention. Kervarrec et al. suggested that VEGFA may be a potential therapeutic target in MCC following promising results from drug trials in mouse models using the monoclonal antibody bevacizumab [38]. The PI3K-Akt and the MAPK signaling pathways are two well-known oncogenic pathways, for which there are numerous established methods of pharmacological inhibition used in the treatment of other malignancies. These include PI3K, Akt, and mTOR inhibitors for the PI3K-Akt pathway and BRAF and MEK inhibitors for the MAPK pathway [39,40]. Furthermore, the insulin-like growth factor pathway constitutes another oncogenic pathway that can be targeted by blocking the IGF-1 receptor or its ligands IGF-1 and IGF-2 [41].
Other notable DAGs included MYB, DST, and KIF23. MYB (c-MYB) encodes an oncogenic transcription factor that regulates processes such as cell proliferation and apoptosis in several other malignancies. Small-molecule inhibitors of c-MYB, such as celastrol and blumbagin, have shown promising results in cell cultures and mouse models [42]. DST encodes a barrier protein that supports melanoma cell growth in vitro and in vivo, likely by interfering with immune cell infiltration or by enhancing angiogenesis [43]. KIF23 encodes a microtubule-associated motor protein involved in the regulation of cytokinesis. Upregulation of KIF23 increases cell proliferation and worsens prognosis in other malignancies, such as gastric cancer. Knockdown of KIF23 resulted in marked inhibition of proliferation in gastric cancer [44]. Other DAGs, the silencing of which has been demonstrated to inhibit cell proliferation or increase apoptosis in other malignancies, include MLF1, MELTF, CIT, RRBP1, CHD7, and MEIS2 [45][46][47][48][49][50].
Another curious finding considering DAGs was enrichment of developmental pathways and processes. This may suggest that the cancer cells of more aggressive MCC revert to a more primitive, stem-cell-like state. An embryonic stem-cell-like gene expression pattern has been found to correlate with poor tumor differentiation and poor prognosis in other malignancies [51]. So-called cancer stem cells that bear specific cell-surface markers and possess the abilities of self-renewal, induction of metastasis, evasion of apoptosis, and resistance to conventional cancer treatments have been described in several other malignancies, but have not as of yet been characterized in MCC. Advances have been made in specifically targeting these cells, including immunotherapy and gene therapy [52].
For SAGs, we found a clear association with pathways and processes related to immune response. Almost all of the most significantly enriched GO BP, GO MF, and KEGG pathway terms were immune-response-related and were especially related to antigen processing, MHC-dependent antigen presentation, and T-cell mediated immune response. These findings underline the importance of a functional immune system, capable of creating a hostile tumor microenvironment, for MCC-specific survival. These findings are also consistent with previous findings on the abundance of tumor-infiltrating lymphocytes, notably CD8-positive T-cells, as a strong predictor of good survival in MCC [6,8,9,11].
The positive correlation between the abundance of CD3-positive, tumor-infiltrating lymphocytes and the expression levels of the SAGs causing the enrichment of the KEGG pathway antigen processing and presentation suggests that the abundance of CD3-positive lymphocytes functions as a surrogate marker for an immunogenic gene expression signature.
Notable SAGs involved in processes not related to immune response included DUSP2, LLGL1, STAT1, and OGFR. DUSP2 encodes a phosphatase that deactivates protumorigenic MAP-kinases, the loss of which predicts a poor prognosis in bladder cancer [53]. LLGL1 encodes a cytoskeleton-associated protein involved in maintaining cell polarity; the loss of LLGL1 is associated with a loss of cellular adhesion, dissemination of cells, and distant metastases in several cancers including gastric cancer and malignant melanoma [54,55]. STAT1 encodes a protein that serves tumor-suppressive functions in many cancers and has been recognized as a potential biomarker for patient selection for treatment with anti-PD-1/anti-PD-L1 antibodies in breast cancer, as p-STAT1 correlates with higher PD-L1 and HLA class I expression [56,57]. Upregulation of opioid growth factor signaling through OGFR (opioid growth factor receptor) suppresses proliferation in several other malignancies, including lung and ovarian cancer [58,59].
Owing to the strong correlation between MCPyV-negativity and MCC-specific death, we repeated the differential expression test of genes based on MCPyV status within the poor prognosis group. Of the 50 DAGs, eight were significantly upregulated in MCPyVnegative tumors, suggesting that their prognostic relevance resulted at least in part from an association with MCPyV-negativity. The remaining 42 DAGs were not significantly differentially expressed based on MCPyV status within the poor prognosis group, suggesting a prognostic relevance regardless of viral status.
It should be noted that the average quality of the extracted RNA was fairly poor (average RIN of 2.1), as is often the case with FFPE samples. However, it has been shown that 3 tag counting, such as that used in this study, markedly decreases the amount of false positives when studying differential expression of genes in samples of varying RIN at the expense of decreased sensitivity [60]. This study utilized a sequencing pipeline optimized for FFPE samples. Specifically, QuantSeq sequences only the 3 end of the RNA transcript, thus significantly reducing the impact of partial RNA fragmentation. This is in contrast to, for example, Illumina's TruSeq, which sequences most of the transcript. In a study comparing Lexogen's QuantSeq and Illumina's TruSeq, there was a strong correlation between the methods concerning the average expression values for all expressed genes. Both methods identified a similar number of expressed protein-coding genes, with QuantSeq identifying approximately 94% of the protein-coding genes found by TruSeq [15]. The correlation analysis between GSE39612 and our transcriptomic data revealed a strong correlation between the logFC-values based on the MCPyV status of the 136 genes studied. This suggests a high specificity of our FFPE data as compared with data from fresh frozen tissues when studying differentially expressed genes. Another challenge with the study design is that, because we analyzed bulk transcriptomic information, we cannot say if a certain expression signature originates from a specific subset of cells in the tumor or if it represents gene expression of the tumor overall.

Conclusions
In summary, we found several novel genes associated with disease-specific survival in MCC. The DAGs were most notably associated with angiogenesis, the PI3K-Akt signaling pathway, the MAPK signaling pathway, and embryonic developmental processes. The SAGs were most notably associated with antigen presentation and immune response.
Further studies are required to determine if some of these genes could be used clinically as prognostic markers, as markers to guide therapy selection, or as targets of precision therapy in MCC.