Identification of Age-Associated Transcriptomic Changes Linked to Immunotherapy Response in Primary Melanoma

Melanoma is a lethal form of skin cancer. Immunotherapeutic agents such as anti-PD-1 (pembrolizumab and nivolumab) and anti-CTLA-4 (ipilimumab) have revolutionized melanoma treatment; however, drug resistance is rapidly acquired. Several studies have reported an increase in melanoma rates in older patients. Thus, the impact of ageing on transcriptional profiles of melanoma and response to immunotherapy is essential to understand. In this study, the bioinformatic analysis of RNA seq data of old and young melanoma patients receiving immunotherapy identifies the significant upregulation of extra-cellular matrix and cellular adhesion genes in young cohorts, while genes involved in cell proliferation, inflammation, non-canonical Wnt signaling and tyrosine kinase receptor ROR2 are significantly upregulated in the old cohort. Several Treg signature genes as well as transcription factors that are associated with dysfunctional T cell tumor infiltration are differentially expressed. The differential expression of several genes involved in oxidative phosphorylation, glycolysis and glutamine metabolism is also observed. Taken together, this study provides novel findings on the impact of ageing on transcriptional changes in melanoma, and novel therapeutic targets for future studies.


Introduction
Malignant melanoma is a lethal form of skin cancer accounting for more than 9000 deaths per year in the US [1]. Despite the implementation of health measures such as the application of sunscreens, the incidence of melanoma is increasing to more than 80,000 cases per year in the US [2]. Melanoma develops due to driver mutations in 48 core melanoma genes [3]. Melanoma development is also driven by mitogen-activated protein kinase (MAPK) pathway activation, which is driven by mutations in BRAF (52%) or NRAS (28%) leading to malignant transformation [3]. The most common BRAF mutations are V600E followed by V600K and V600R [3]. Resistance to BRAFV600E inhibitors after treatment remains a challenge to melanoma treatment. At later stages, melanoma resection remains the main treatment with poor prognosis [4]. The development of immunotherapy has revolutionized melanoma treatment; however, melanoma acquires resistance, making treatment more challenging [5]. Aging is an inevitable biological process and a risk factor for many diseases, including melanoma. Aging is associated with the accumulation of damaged molecules leading to replicative senescence in dividing cells, telomere shortening and reduced levels or the absence of telomerase (hTERT) expression [6]. Senescent cells constantly produce proinflammatory cytokineschemokines-growth factors leading to a microenvironment favorable for the development of age-related diseases including cancer. Several studies have shown that the incidence of melanoma increases with age. Melanoma rates increased by an annual percentage change of 1.8% in adults aged 40 and higher [7]. In addition, there has been a statistically significant decrease in melanoma rates in young adults [7]. Patients aged 65 and older are the most vulnerable age group affected by melanoma-related deaths [8]. Animal experiments have shown that the melanoma cell line Yumm1.7 derived from a BrafV600E/Cdkn2a−/−/Pten−/− mouse model of melanoma grows slowly but aggressively in aged mice with increased angiogenesis and metastases due to the age-related increase in the Wnt antagonist sFRP2 [9]. Luckily, old melanoma patients have been found to have a better response to the immune checkpoint inhibitor anti-PD1 due to accumulation of a FOXP3+ Treg cell population [10]. It has also been observed that in young mice as well as young adults with melanoma, there is an accumulation of a significantly higher population of Treg cells driving resistance to immune check point inhibitors [10]. It is, therefore, important to comprehensively understand the molecular mechanisms driving melanoma progression during aging to identify new therapeutic targets for the treatment of malignant melanoma. Elucidating the transcriptional profiles of aging and their association with physiological and pathological changes in melanoma is, therefore, critical to identify new therapeutic targets. High-throughput genome-wide transcriptomic profiling is a powerful tool that identifies gene expression signatures. The Cancer Genome Atlas (TCGA) is a comprehensive project cataloguing genomic and transcriptomic data of cancers by the National Cancer Institute (NCI), the National Human Genome Research Institute (NHGRI) and the National Institute of Health (NIH) [11]. Since transcriptomic states of aged and young melanoma patients are not fully understood, leading to a knowledge gap, in this study, robust bioinformatic analysis of publicly available RNA seq data from 13 patients receiving immunotherapy (representative of two distinct age groups) was performed. Five old and eight young patients diagnosed with melanoma were collected and an analysis of gene expression profiles was performed. Transcriptomic changes and gene expression signatures associated with response to immunotherapy were investigated.

Clinical Cohorts
Two cohorts of old (n = 5) and young (n = 8) immunotherapy-treated patients (ipilimumab) diagnosed with primary malignant melanoma were identified from TCGA and served as datasets for the evaluation of age-associated transcriptomic changes. The class of phenotypes that was used was primary tumor in the clinical category of "primary tumor".

Data Selection and Processing
Transcriptomic RNAseq gene expression level 3 data containing reads per kilobase per million mapped reads (RPKM) were downloaded for each case. RPKM is a widely used RNAseq normalization method and is computed as follows: RPKM = 109(C/NL), where C is the number of reads mapped to the gene, N is the total number of reads mapped to all genes, and L is the length of the gene. RPKM results were generated with SeqWare pipeline. The reference genome was genome GRCh38.p0, and the genome name is GRCh38.d1.vd1. The alignment of raw reads was completed with BWA.

Differential Expression Analysis of Individual Genes
Differentially expressed genes (DEGs) were analyzed using the DESeq2 Bioconductor R package [12]. Normalization was based on the relative log expression method in DESeq2. The identification of significant DEGs was based on a log 2 fold change cutoff value of ≤2 (for downregulated genes) or ≥2 (for upregulated genes). A significance level of an adjusted p value of <0.05, using a false discovery rate (FDR) cutoff of <0.1, was set. Heatmaps were generated using "pheatmap", "RColorBrewer", "ComplexHeatmap" and "circlize" R packages [13]. A volcano plot was visualized using the EnhancedVolcano R package [14].

Functional and Pathway Enrichment Analysis
Entrez-ID of each DEG was obtained with the "org.Hs.eg.db" R package, then Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the "clusterProfiler," "enrichplot," and "ggplot2" R packages. KEGG analysis is used to discover pathways enriched in genes in a gene set of interest. Gene set enrichment analysis (GSEA), which is not restricted by DEGs, was also used.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Mapping Analysis
The Database for Annotation, Visualization and Integrated Discovery (DAVID) was used to analyze the differentially expressed genes. Pathways were also analyzed by the KEGG program. In KEGG pathway enrichment analysis, enriched pathways were identified according to p < 0.05.

Principal Component Analysis
Principal component analysis was performed with the R packages FactoMineR and Factoextra. Two-dimensional PCA plots were generated using the R ggplot2.

Statistical Analyses
The significance level of the adjusted p value was set at <0.05 and a false discovery rate (FDR) cutoff of <0.1 was set. All analyses were performed with R version 4.1.1 (10 August 2021) using the following packages: gplots, RColorBrewer, org.Hs.eg.db, statmod, DESeq2, EnhancedVolcano, genefilter, pheatmap, NMF, ggplot2, scales, viridis, gridExtra, reshape2, ggdendro, IHW and dendextend. All graphs were generated in R version 4.1.1, R Core Team (2022), R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria.

Demographic and Clinical Features of the Patients
Patients with primary diagnosis malignant melanoma (disease type: nevi and melanomas) were classified into two groups according to age. The young age group had an average age at diagnosis of 42.8 ± 7.75 years. The old group had an average age at diagnosis of 70 ± 4.5 years, p = 0.00001 ( Figure 1A). None of the patients displayed prior or synchronous malignancies. Survival analysis showed that the old age group had significantly worse melanoma-specific survival, p = 0.006 ( Figure 1B).

Identification of Age-Associated Differentially Expressed Genes in Melanoma
The quality of RNA-seq data was evaluated with principal component analysis (PCA) and the clustering of RNA-seq samples using Euclidean distance. PCA analysis plots of the 13 RNAseq samples showed a clear separation of samples along principal component 1 (PC1), which explained around 70% of total variance ( Figure 1C). In addition, Euclidean sample distances showed that the two age groups were well-clustered ( Figure 1D). Differential expression analysis was performed using the DESeq2 R package to identify upregulated and downregulated genes between the old and young patients diagnosed with primary malignant melanoma. DEGs were visualized as an MA plot (log 2 fold change vs. mean of normalized counts) of young vs. old. Red dots represent transcripts with positive and negative log 2 fold change values, and also indicate the upregulation and downregulation of DEGs (Figure 2A). A total of 3112 significant DEG genes were identified using a statistical cutoff: q < 0.05. A total of 1345 genes were upregulated, and 1767 genes were downregulated in the young cohort (Figure 2A,B). The analysis of DEGs using a hierarchical clustering heatmap revealed distinct gene expression profiles between the two age groups ( Figure 3A). At cutoff log 2 fold change values of 3 and −3, there were 200 upregulated genes and 386 downregulated genes in the young cohort compared to the old cohort, respectively. The genes with the lowest p values were ATPase family AAA domain-containing 3B (ATAD3B) and immunoglobulin-like and fibronectin type III domaincontaining 1 (IGFN1). Microsomal glutathione S-transferase 1 (MGST1), a gene that plays a role in inflammation, was remarkably upregulated in the young age group (log 2 fold change = 3.02) ( Figure 3B). Four collagen genes involved in the extracellular matrix and cellular adhesion were upregulated in the young group: COL14A1, COL2A1, COL6A6, COL4A4 ( Figure 3B). We also observed a significant upregulation of cell proliferation genes p21CIP1 (CDKN1A) and CDK2 in the young cohort, while CDK1 was upregulated in the old cohort p < 0.005 ( Figure 3C). A slight upregulation of ARMC8 was also observed in the young cohort ( Figure 3C). ARMC8 plays critical roles in cell proliferation, apoptosis, and differentiation, and has been reported to be a prognostic marker in several malignancies including liver, lung and breast cancers [15][16][17][18]. It has been reported that ARMC8 is upregulated in malignant melanoma cell lines and is associated with an increased invasiveness of melanoma [15]. We also observed the significant upregulation of interleukin-17A (IL-17A) and interleukin 11 (IL11) in the old cohort p < 0.05 ( Figure 3D). Interleukin-17A (IL-17A) is produced by Th17 cells which infiltrate the tumor microenvironment and induces the expression of several inflammatory cytokines including IL-1β, IL-16 and IL-23, leading to variable effects on tumor growth at different stages [19,20]. IL-17A upregulation has been reported in a 53-year-old male with collision primary laryngeal malignant melanoma and invasive squamous cell carcinoma [21]. IL-11 is a pleiotropic interleukin involved in tumor development. IL-11 has been shown to be elevated in colorectal cancer and in the exosomes of metastatic uveal melanoma and is associated with poor prognosis in melanoma [22][23][24][25]. It is known that phenotype switching is characteristic of melanoma, where proliferative cells tend to be non-invasive, while invasive cells tend to be non-proliferative, a process regulated by the Wnt signaling pathway [26][27][28]. We identified the dysregulation of several genes in the Wnt signaling pathway ( Figure 3E). Interestingly, the non-canonical Wnt molecule Wnt5A and the tyrosine kinase receptor ROR2 which drive invasion, metastasis and therapeutic resistance [29][30][31], were significantly upregulated in the old cohort ( Figure 3E).

Identification of Age-Associated Differentially Expressed Genes in Melanoma
The quality of RNA-seq data was evaluated with principal component analysis (PCA) and the clustering of RNA-seq samples using Euclidean distance. PCA analysis plots of the 13 RNA-seq samples showed a clear separation of samples along principal component 1 (PC1), which explained around 70% of total variance ( Figure 1C). In addition, Euclidean sample distances showed that the two age groups were well-clustered ( Figure 1D). Differential expression analysis was performed using the DESeq2 R package to identify upregulated and downregulated genes between the old and young patients diagnosed with primary malignant melanoma. DEGs were visualized as an MA plot (log2 fold change vs. mean of normalized counts) of young vs. old. Red dots represent tran- vated in colorectal cancer and in the exosomes of metastatic uveal melanoma and is associated with poor prognosis in melanoma [22][23][24][25]. It is known that phenotype switching is characteristic of melanoma, where proliferative cells tend to be non-invasive, while invasive cells tend to be non-proliferative, a process regulated by the Wnt signaling pathway [26][27][28]. We identified the dysregulation of several genes in the Wnt signaling pathway ( Figure 3E). Interestingly, the non-canonical Wnt molecule Wnt5A and the tyrosine kinase receptor ROR2 which drive invasion, metastasis and therapeutic resistance [29][30][31], were significantly upregulated in the old cohort ( Figure 3E).

Pathway Analysis, Functional Annotation of Differentially Expressed Genes and Identification of Key Regulatory Genes of Response to Immunotherapy
Functional annotation clustering of the most upregulated and downregulated genes (log2 fold change = 3 and −3, respectively) was performed with DAVID to identify the biological processes that are related to differential gene expression changes in the two melanoma age groups receiving immunotherapy. Among the 200 genes upregulated in the young group, the top clusters were related to signal peptide, glycoprotein, glycosylation site:N-linked (GlcNAc) and ion transport. The top clusters in the old age group were related to extracellular region, signal peptide, extracellular space and glycosylation site:N-

Pathway Analysis, Functional Annotation of Differentially Expressed Genes and Identification of Key Regulatory Genes of Response to Immunotherapy
Functional annotation clustering of the most upregulated and downregulated genes (log 2 fold change = 3 and −3, respectively) was performed with DAVID to identify the biological processes that are related to differential gene expression changes in the two melanoma age groups receiving immunotherapy. Among the 200 genes upregulated in the young group, the top clusters were related to signal peptide, glycoprotein, glycosylation site:N-linked (GlcNAc) and ion transport. The top clusters in the old age group were related to extracellular region, signal peptide, extracellular space and glycosylation site:N-linked (GlcNAc). Pathway and gene set enrichment analysis of the upregulated genes showed the significant enrichment of genes involved in several biological processes, most notably, interferon-α, interferon-γ response, IL2-STAT5 signaling, JAK-STAT3, P53 and inflammatory response pathways ( Figure 4A). KEGG pathway analysis of DEG genes was also performed, and it was revealed that the pathways that had the highest enrichment scores for the DEG genes were riboflavin metabolism (0.71), histidine metabolism (0.64), antigen processing and presentation (0.6), Toll-like receptor signaling, (0.6), circadian rhythm (0.57), lysosome (0.55), calcium signaling and the PD-L1 expression/PD-1 checkpoint pathway in cancer ( Figure 4B,C). It is known that cancer cells upregulate programmed death ligand-1 (PD-L1) which binds inhibitory receptor programmed death receptor-1 (PD-1) on the T cell surface to avoid immune attack [32]. Several studies have evaluated the efficacy of PD-1 and PD-L1 inhibitors in different cancers, including lung cancer, renal cancer and malignant melanoma in PD-L1-negative and PD-L1-positive tumors. Therapeutic benefit from immunotherapy was observed in patients with the PD-L1-negative tumors; however, resistance to immunotherapy remains a challenge [33,34]. Indeed, gene set enrichment analysis (GSEA) of transcriptomic data confirmed the significant enrichment of genes in the T cell receptor signaling, PD-L1 expression and PD-1 checkpoint and JAK-STAT signaling pathways, which are essential pathways in melanoma ( Figure 4D-G) [35][36][37].

Differential Expression of Treg Signature Genes and Metabolic Genes in Melanoma
FOXP3+ T regulatory cells (Tregs) are characterized by the expression of the transcription factor forkhead box P3 (FOXP3). They infiltrate the tumor microenvironment, thus providing an immunosuppressive microenvironment and immune tolerance [38]. Functionally immunosuppressive Tregs are also highly enriched in the tumor microenvironment of patients with melanoma [39][40][41]. Tregs provide an immune-tolerant microenvironment by secreting a number of factors such as IL-10, TGF-β, IL-35, PD-1, LAG-3 and TIM-3 [42,43]. The inhibition of Treg function is therefore required for tumor clearance and therapeutic effectiveness. It has been shown in a cohort of primary melanoma samples from NYU and Vanderbilt that young patients acquire resistance to immunotherapy due to a reduction in FOXP3+ Treg cell population, while old patients (> 66 years) have a small CD8+ T cell population [10]. Animal experiments in murine models of melanoma have shown that the depletion of Tregs leads to tumor clearance and increased survival, indicating the critical importance of this cell type for antitumor immunity [44]. Thus, in this analysis, Treg signature genes in young and old patients treated with the immunotherapeutic agent ipilimumab were compared to dissect the molecular pathways regulating T cell function in melanoma during aging. Two Treg signature genes, NR4A3 and IKZF2, were significantly upregulated (log 2 fold change > 1.5) ( Figure 5A). Nr4a3 is an orphan nuclear receptor upregulated in T cells undergoing differentiation and is required for FOXP3 induction in T cells [45][46][47][48]. IKZF2 (Helios) is a transcription factor required for the suppression of IL-2 production in Tregs and regulates FOXP3 binding to the IL2 promoter, thus silencing Il2 transcription in Tregs. [49]. The expressions of the most prominent signature genes of tumor-infiltrating Tregs, such as CXCR5, IL17F, IL17A, IL22, FOXP3, IL12RB2, TNFRSF9, CD274, TNFRSF4, TNFRSF9, IL10, CCR7 and STAT1, were also investigated ( Figure 5A). ID3, RDH10 and EOMES are transcription factors associated with a dysfunctional tumor-infiltrating T cell state [50].
ID3 and RDH10 were enriched in the old cohort, while EOMES was enriched in the young cohort ( Figure 5A). We also observed a significant upregulation of STAT1 in the young cohort ( Figure 5A). STAT1 is involved in cytokine production and is required for the recruitment and activation of T cells in the tumor microenvironment. STAT1 has also been found through a two-cell CRISPR-type screen including human T cells as effectors and melanoma cells as targets to be upregulated, affecting immunotherapy response [51]. The GTP-binding protein 4 (GBP4) was previously identified as a prognostic signature gene that separates melanoma patients into low-and high-risk groups according to survival [52]. GBP4 was significantly upregulated in the young group ( Figure 5A). In addition, we examined metabolic gene expression in the young and aged cohorts. It is well known that metabolic rewiring is a hallmark of all types of cancer [53][54][55]. It has been previously shown that immunotherapy (PD-1 blockade)-resistant melanoma tumor cells acquire a hypermetabolic state by the upregulation of glycolysis and mitochondrial oxidative phosphorylation. The upregulation of lactate and TCA cycle metabolites is observed in immunotherapy-resistant melanoma cells [56]. Glutamine addiction is a hallmark of malignant transformation, including melanoma [57]. Glutamine is an anaplerotic amino acid. Carbon derived from glutamine is used to maintain TCA cycle intermediates and its nitrogen is used for transamination reactions, as well as purine and redox intermediates (NAD and NADP) biosynthesis [58]. The upregulation of genes involved in the biosynthesis of proline from glutamate is also a characteristic of malignant melanoma [59,60] and the depletion of glutamine sensitizes melanoma cells to TRAIL-mediated cell death [61]. Five genes-ATP12A, ATP6AP1, ATP6V0A1, ATP6V0C and ATP6V0D2-involved in the oxidative phosphorylation pathway were upregulated in the young cohort ( Figure 5B), while six genes-ATP6V1C2, COX6C, COX7B, PPA1, UQCRB and UQCRH-were upregulated in the aged cohort ( Figure 5B). Glycolysis genes ALDH9A1, LDHAL6A, PFKL and PKLR were upregulated in the young cohort ( Figure 5B), while the glutamine metabolism gene GLS2 was notably upregulated in the aged cohort ( Figure 5B). The dysregulation of (NFE2L2 or NRF2) target genes was also observed. The transcription factor nuclear factor erythroid 2-related factor 2 (NFE2L2 or NRF2) is the master regulator of cellular redox homeostasis. NRF2 regulates the expression of genes involved in antioxidant defense and xenobiotic metabolism. NRF2 is activated in melanoma and is required for melanoma cell proliferation; its knockdown sensitizes melanoma cells to oxidative stress [62]. NRF2 expression was slightly upregulated in the young cohort while its target genes NQO1 and ALAS1 were significantly upregulated in the young cohort. GSTT2B, GSTT2 and TXNDC17 were significantly upregulated in the old cohort ( Figure 5B). Finally, an analysis of the correlation between gene expression and melanoma patients' survival from TCGA showed that CCR7 (hazard ratio (HR) 0.67, p = 0.0029), CXCR5 (HR 0.7, p = 0.01), EOMES (HR 0.76, p = 0.047), GBP4 (HR 0.5, p = 4.9 × 10 −7 ), TNFRSF9 (HR 0.59, p = 0.00011) and STAT1 (HR 0.58, p = 5.3 × 10 −5 ) were significantly related to patient prognosis. The overall survival rate of patients with high expressions of CCR7, CXCR5, EOMES, GBP4, TNFRSF9 and STAT1 in melanoma was significantly higher ( Figure 5C-H), which is consistent with the significantly higher survival rate observed for the young cohort receiving immunotherapy ( Figure 1B).

Discussion
Melanoma accounts for 80% of skin cancer deaths [63]. Although immunotherapy has revolutionized melanoma treatment, many melanoma patients remain resistant to immunotherapy [64,65]. Immune checkpoint inhibitor proteins such as CTLA-4 and PD-1 are expressed in activated T cells. CTLA-4 inhibits the activation of T cells and PD-1 binds to PD-L1/PD-L2 ligands that are expressed in melanoma, thus inhibiting immune attack [63]. Aging is a risk factor for melanoma incidence. In this study, the impact of aging on the transcriptomic profiles of melanoma patients receiving immunotherapy was investigated. A comprehensive analysis, including gene expression, the identification of DE genes, GSEA and pathway enrichment, was performed to identify signature genes. Our

Discussion
Melanoma accounts for 80% of skin cancer deaths [63]. Although immunotherapy has revolutionized melanoma treatment, many melanoma patients remain resistant to immunotherapy [64,65]. Immune checkpoint inhibitor proteins such as CTLA-4 and PD-1 are expressed in activated T cells. CTLA-4 inhibits the activation of T cells and PD-1 binds to PD-L1/PD-L2 ligands that are expressed in melanoma, thus inhibiting immune attack [63]. Aging is a risk factor for melanoma incidence. In this study, the impact of aging on the transcriptomic profiles of melanoma patients receiving immunotherapy was investigated. A comprehensive analysis, including gene expression, the identification of DE genes, GSEA and pathway enrichment, was performed to identify signature genes. Our results show that the old-age cohort had significantly worse melanoma-specific survival p < 0.001 ( Figure 1A). A total of 1345 DE genes were upregulated and 1767 DE genes were downregulated in the young cohort (Figure 2A,B). The results reveal signatures of the impact of age on melanoma. MGST1 plays a critical role in inflammation, is overexpressed in cancer, and correlates with drug resistance [66]. In our analysis, MGST1 was remarkably upregulated in the young cohort. Prall et al. reported that MGST1 expression is age-dependent and Zeng et al. reported the overexpression of MGST1 in high-risk melanoma patients [67,68]. Collagen genes COL14A1, COL2A1, COL6A6, COL4A4 involved in extracellular matrix organization were upregulated in the young group ( Figure 3B). This is consistent with previous studies that have reported significant changes in the skin's extracellular matrix during aging, such as collagen loss [69]. Miskolczi et al. also reported that melanoma cell adhesion and nuclear YAP localization are regulated by the mechanical properties of collagen. Collagen stiffness induced the expression of melanoma differentiation genes TRPM1, PMEL, TYR and MLANA, as well as well as the proliferation and survival genes CDK2 and BCL2A1 [70]. Genes regulating cell proliferation, apoptosis and differentiation were upregulated in the young cohort ( Figure 3C). The activation of inflammatory genes such as IL-17A and IL-11 was observed in the old cohort ( Figure 3D). IL-17A induces the expression of inflammatory cytokines IL-1β, IL-16 and IL-23 in the tumor microenvironment in malignant melanoma [19][20][21]. The activation of such inflammatory pathways is potentially associated with poor survival observed for the old cohort receiving immunotherapy. Indeed, Mehta et al. have reported that the inflammation-induced dedifferentiation of melanoma cells contributes to poor survival and resistance to immunotherapy in aged melanoma patients [71]. It has also been reported that phenotype switching in melanoma is regulated by the Wnt signaling pathway and that Wnt5A-treated B16 melanoma cells acquire greater metastatic and invasive potential, an effect mediated by the orphan tyrosine kinase receptor ROR2 [30,72]. Indeed, Wnt5A and ROR2 were upregulated in the old cohort receiving immunotherapy ( Figure 3E). This is consistent with previous results indicating that Wnt5A expression increases in old (55-65 years old) melanoma patients compared to young (25-35 years old) patients [73]. Behera et al. also reported that melanoma resistance to vemurafenib is driven by Wnt5A and that PPARγ activation with rosiglitazone upregulated the age-related protein klotho, a fibroblast growth factor-23 (FGF23) receptor, and decreased Wnt5A expression in therapy-resistant old melanoma patients, thus reducing tumor burden [74]. Our results further reveal the enrichment of genes in several biological processes, including interferon-α, IL2-STAT5 signaling, and the JAK-STAT3 and P53 pathways ( Figure 4A). Interestingly, KEGG pathway analysis of DE genes showed the upregulation of PD-L1 expression/the PD-1 checkpoint pathway in cancer ( Figure 4B-D). DE gene analysis also revealed the differential expression of several Treg signature genes and metabolic genes. Immunosuppressive FOXP3+ Tregs infiltrate the tumor microenvironment, inhibiting immune attack and thus affecting the response to immunotherapy [10,75,76]. Treg signature genes CXCR5, TNFRSF9, CCR7, STAT1, GBP4, and EOMES were significantly upregulated in the young cohort and were correlated with higher survival in melanoma ( Figure 5B-H), while ID3, IL-17A, IL-17F, RDH10 and IL-11 were significantly upregulated in the old cohort ( Figure 5A). Indeed, aging impacts the infiltration of immune cells and the expression of immune checkpoints leading to an immunosuppressive microenvironment and an aggressive phenotype of tumor cells [77]. We also observed the differential expression of several genes involved in metabolic rewiring in melanoma. Oxidative phosphorylation genes ATP12A, ATP6AP1, ATP6V0A1, ATP6V0C and ATP6V0D2 were upregulated in the young cohort ( Figure 5B), while ATP6V1C2, COX6C, COX7B, PPA1, UQCRB and UQCRH were upregulated in the aged cohort ( Figure 5B). GLS2 was upregulated in the aged cohort. It has been reported that GLS2 upregulation increases tumor metastasis and correlates with poor survival [78]. Consistent with our transcriptomic analysis, Wu et al. reported that aging strongly impacts several biological pathways in cancer, including epithelial-mesenchymal transition (EMT), metabolism, KRAS signaling, inflammatory response, glycolysis and Il2-stat5 signaling [77].

Conclusions
In conclusion, these data reveal transcriptional changes associated with age; the dysregulation of many biological, inflammatory and metabolic pathways associated with greater metastatic and invasive potential in old melanoma patients receiving immunotherapy; and potential novel therapeutic targets in melanoma. One limitation in the current study is the small sample size due to limited data on melanoma patients receiving immunotherapy in the TCGA database. In the future, large samples, clinical cohorts and basic experiments will be required to verify the impact of ageing on immunotherapy response and the prognosis of melanoma patients.

Data Availability Statement:
The datasets of the current study analyzed with R are available from the corresponding author on request.

Conflicts of Interest:
The author declares no conflict of interest.