Identification of 15 lncRNAs Signature for Predicting Survival Benefit of Advanced Melanoma Patients Treated with Anti-PD-1 Monotherapy

The blockade of programmed cell death protein 1 (PD-1) as monotherapy has been widely used in melanoma, but to identify melanoma patients with survival benefit from anti-PD-1 monotherapy is still a big challenge. There is an urgent need for prognostic signatures improving the prediction of immunotherapy responses of these patients. We analyzed transcriptomic data of pre-treatment tumor biopsies and clinical profiles in advanced melanoma patients receiving only anti-PD-1 monotherapy (nivolumab or pembrolizumab) from the PRJNA356761 and PRJEB23709 data sets as the training and validation cohort, respectively. Weighted gene co-expression network analysis was used to identify the key module, then least absolute shrinkage and selection operator was conducted to determine prognostic-related long noncoding RNAs (lncRNAs). Subsequently, the differentially expressed genes between different clusters were identified, and their function and pathway annotation were performed. In this investigation, 92 melanoma patients with complete survival information (51 from training cohort and 41 from validation cohort) were included in our analyses. We initiallyidentified the key module (skyblue) by weighted gene co-expression network analysis, and then identified a 15 predictive lncRNAs (AC010904.2, LINC01126, AC012360.1, AC024933.1, AL442128.2, AC022211.4, AC022211.2, AC127496.5, NARF-AS1, AP000919.3, AP005329.2, AC023983.1, AC023983.2, AC139100.1, and AC012615.4) signature in melanoma patients treated with anti-PD-1 monotherapy by least absolute shrinkage and selection operator in the training cohort. These results were then validated in the validation cohort. Finally, enrichment analysis showed that the functions of differentially expressed genes between two consensus clusters were mainly related to the immune process and treatment. In summary, the 15 lncRNAs signature is a novel effective predictor for prognosis in advanced melanoma patients treated with anti-PD-1 monotherapy.


Introduction
Melanoma is one of the most aggressive malignant skin tumors, and the incidence has been increasing worldwide in recent decades [1,2]. Since the US Food and Drug Administration (FDA) and European Medicines Agency (EMA) approved a variety of immune checkpoint inhibitors (ICI) therapies for advanced melanoma in 2011, the overall mortality of advanced melanoma fell by 17.9% from 2013 to 2016 [3], but it is still at a very high level [4].
During the last years, biomarkers for effectiveness of tumor immunotherapy, including genomic instability as described by microsatellite instable (MSI) or tumor mutational burden (TMB) status, and immune cell infiltration into tumors, have been defined [5]. As one of the cancer types with high TMB and high immune cell infiltration, melanoma can be considered for immunotherapy [5]. Despite all efforts of early diagnosis, metastatic melanoma still has a poor prognosis and remains a challenge for treating physicians. Existing ICI therapies include the blockade of programmed cell death protein 1 (PD-1)/programmed deathligand 1 (PD-L1) and the cytotoxic T-lymphocyte antigen 4 (CTLA-4) pathway. Therapies based on blockade of PD-1 in human melanoma achieved a success [6,7]; nivolumab (an anti-PD-1 monoclonal antibody) [8], and pembrolizumab (another anti-PD-1 monoclonal antibody) [9] monotherapy improve relapse/recurrence-free survival of stage III melanoma patients. Moreover, anti-PD-1 monoclonal antibodies indicate superior overall/recurrencefree survival versus anti-CTLA-4 agent for advanced melanoma [8,10,11]. Melanoma patients treated with anti-PD-1 monotherapy have a longer progression-free survival (PFS) compared to those treated with both strategies [12]. Further, anti-PD-1 monotherapy provides survival benefits in responding patients, but still many patients fail to respond [12]. Identifying the responsive population is a priority to optimize drug selection and improve patient outcomes. Future research should focus on identifying additional biomarkers to select patients who are most likely to benefit from certain immunotherapies.
Numerous long noncoding RNAs (lncRNAs) have been identified in human genomes [13] and they are nowadays considered as prognosis signatures in various tumors [14,15]. However, there has been limited systematic characterization of these elements in melanoma, especially in melanoma patients treated with anti-PD-1 monotherapy. Here, we identified a 15 lncRNAs signature in advanced melanoma by weighted gene co-expression network analysis (WGCNA) and logistic least absolute shrinkage and selection operator (LASSO) in the training cohort (PRJNA356761) and validated it in the validation cohort (PRJEB23709) to determine and predict the response of melanoma patients receiving anti-PD-1 monotherapy.

Module Construction
After crossing all lncRNAs in the two cohorts, unsigned co-expression modules were constructed in the training cohort using the WGCNA algorithms in R as described previously [20]. We used one-step network construction method to identify co-expression modules through the blockwiseModules function in the WGCNA package [21]. Before coexpression network construction, the flashClust tool in R was utilized to perform hierarchical clustering analysis of the samples with the appropriate threshold value to detect and eliminate the outliers. According to scale-free topology criterion, a soft-thresholding power β (the power values ranged from 1 to 20), which was calculated by the pickSoftThreshold function of the WGCNA, was chosen to build an adjacency matrix [21]. In our study, the power of 6 was used for this network.
Then, the topological overlap matrix was constructed based on the adjacency matrix. A dissimilarity matrix was used to detect gene modules (gene sets with high topological overlap) through a dynamic tree cutting algorithm [22,23]. To obtain moderately sized modules, the minimum number of genes was set at 30 and a cutline was chosen to merge modules with similar expression patterns. To identify the relationships between modules and clinical traits, we calculated the correlation between module eigengenes and clinical trait and searched for the most significant associations. The module eigengenes was calculated by the first principal component, which were considered as a representative of the expression patterns of module genes [24]. For each module, we defined the module membership as the correlation of gene expression profile with module eigengene and the gene significance as the absolute value of the correlation between gene and clinical traits. In this study, genes with high module membership in a module were assigned to the module and the module with high gene significance and p value < 0.05 was considered to be highly related to clinical traits. Moreover, after introducing validation cohort, we used network-based statistics, which generated a composite statistic value (Z summary ) using a permutation test to measure the strength of lncRNAs module and expression module preservation, to assess whether the density and connectivity patterns of lncRNAs were also preserved [25]. Z summary > 5 implies strong evidence for module preservation [26,27]. Since the Z summary statistic bias towards a module with a large size [25], a rank-based statistic medianRank, calculated from observed preservation values and conducted no permutation test against background gene modules, was used to measure the relative preservation irrespective of module size [28].

Identifcation of lncRNAs Signature
Although which module was most relevant to the predictor of melanoma can be identified after WGCNA, we further applied a novel, modern statistical shrinkage technique to examine the association between lncRNAs and the prognosis of melanoma to establish prognostic lncRNAs signature. The logistic LASSO regression model is a shrinkage method that can actively select from a large and potentially multicollinear set of variables in the regression, resulting in a more relevant and interpretable set of predictors [29]. One interesting property of LASSO is that the estimates of the regression coefficients are sparse, which means that many components are exactly 0 [30]. We utilized the glmnet package (version 2.0-16) to fit the logistic LASSO regression.

Development and Validation of the lncRNAs Signature
After the identification of the predictive module, we further clustered melanoma population into different consensus clusters though an optimum cutoff value identified by survivalROC package (version 1.0.3). Ultimately, the lncRNAs signature can distinguish different consensus clusters, and the survival package was applied to perform Kaplan-Meier analysis with the log-rank test to analyze the overall survival (OS) and PFS in both training and validation cohorts. Moreover, we used a time-dependent receiver operating characteristic (ROC) curve to assess the survival prediction, and the area under the ROC curve (AUC) value was computed with the ROCR package (version 1.0.-7) to measure prognostic or predictive accuracy, as described previously [31]. In addition, we calculated the response rates of anti-PD-1 therapy based on the lncRNAs signature in both, training and validation cohort.

Functional Analysis
Subsequently, the differentially expressed genes (DEGs) between different consensus clusters distinguished by the lncRNAs signature were identified by limma package [32] (version 3.28.14) with the cutoff of p < 0. Functional analysis were performed using the clusterProfiler package [33] (version 3.18.0) and msigdbr package [34] (version 7.2.1) to expand our understanding of those lncRNAs signature-related functions and their coordinated regulatory networks.

Immune Cell Enrichment Analysis by xCell
Immune cell enrichment analysis was conducted by xCell function [35] in immunedeconv package (version 2.0.3) [36].

Patient Characteristics
In total, 100 tumor biopsies from melanoma patients treated with anti-PD-1 monotherapy (nivolumab) were assessed by whole-exome sequencing in PRJNA356761, among them, 49 biopsies were pre-therapy [16]. PRJEB23709 included 158 tumor biopsies from 120 melanoma patients treated with anti-PD-1 monotherapy (nivolumab or pembrolizumab) or combined anti-PD-1 and anti-CTLA-4 (nivolumab or pembrolizumab combined with ipilimumab) [12]. In our study, we included 51 melanoma patients form PRJNA356761 and 41 melanoma patients form PRJEB23709 (Figure 1). The patient characteristics are listed in Table 1. Responders were defined as patients with a Response Evaluation Criteria in Solid Tumors (RECIST) response [37] of complete response (CR), partial response (PR), or stable disease (SD) of greater than 6 months with no progression, and non-responders as progressive disease (PD) or SD for less than or equal to 6 months before disease progression.

WGCNA and Key Module Identification
We both obtained data of 16,899 lncRNAs from the training and validation cohorts, by transcriptome analysis, and 4653 lncRNAs after the intersection were included to construct co-expression networks by WGCNA. To exclude the outliers, we chose 85 for the cut tree height for the samples (Figure 2A), and 48 samples were included for subsequent analysis. Then we identified the soft-thresholding power β of six to construct a scalefree network ( Figure 2B,C). As a result, 12 co-expression modules were identified by gathering similarly expressed lncRNAs ( Figure 3A). Interactions between 12 modules were subsequently analyzed. The heatmap demonstrated the topological overlap matrix among all of 4653 lncRNAs in our study ( Figure 3B), indicating that each module showed independent validation to each other. Then, the correlations between module eigengene, and clinical traits were discovered. Moreover, after introducing the validation cohort, we plotted the preservation median rank and Z summary for the modules as a function of module size. The 12 modules (gold, yellow, pink, turquoise, magenta, green, lightgreen, white, darkred, skyblue, skyblue3, and lightsteelblue1) showed strong evidence of preservation (Z summary > 5) ( Figure S1). The module eigengenes in the skyblue module showed a higher correlation with PFS status, OS status, OS, PFS, and responder (R PFS status 2 = −0.43, p < 0.002; R OS status 2 = −0.27, p < 0.05; R OS 2 = 0.24, p < 0.1; R PFS 2 = 0.48, p < 0.0001; R responder 2 = 0.29, p = 0.04; respectively) ( Figure 3C). We therefore chose the skyblue module for further analyses.

Identification of lncRNAs Signature
There were 63 lncRNAs in the skyblue prognosis module we identified (Table S1) (Table 2).

Development and Validation of the lncRNAs Signature
On the basis of the time-dependent ROC curve analysis, the optimal cutoff value that could be used for the 15 lncRNAs signature to stratify melanoma patients treated with anti-PD-1 monotherapy into the high-or low-risk group was determined to be 0.39 in the training cohort ( Figure 4A-C). All 51 patients in the training cohort were segregated into the high-risk group (n = 30) and the low-risk group (n = 21), and the low-risk group exhibited significantly better OS than the high-risk group (hazard ratio (HR) = 0.2855, 95% confidence interval (CI) = 0.1302~0.6259, p = 0.0010, Figure 5A). For the low-risk group, the median OS was not reached, whereas the median OS was 12.5 months (95% CI = 7.46~21.2) for the high-risk group ( Figure 5A). Regarding PFS, all 51 patients in the training cohort were similarly segregated into the low-and high-risk groups, and the low-risk group was correlated with significantly favorable PFS (HR = 0.2805, 95% CI = 0.1427~0.5514, p = 0.0001, Figure 5B). For the low-risk group, the median PFS was 10.95 months (95% CI = 7.20~NA), whereas the median PFS was 1.87 months (95% CI = 1.71~5.36) for the high-risk group ( Figure 5B).    To examine the robust and realistic application of the lncRNAs signature, the performance of the 15 lncRNAs signature was validated in the validation cohort ( Figure 4D-F). The developed 15 lncRNAs signature could actively predict OS and PFS in melanoma patients treated with anti-PD-1 monotherapy in the validation cohort. The lncRNAs signature significantly stratified patients into low-and high-risk groups in terms of OS; more specifically, all 41 patients were segregated into the low-risk group (n = 22) and the high-risk group (n = 19) and showed significantly different OS rates (HR = 0.4634, 95% CI = 0.2054~1.045, p = 0.0580, Figure 5C) according to the optimum cutoff point (−0.321) acquired from the training cohort ( Figure 4D). For the low-risk group, the median OS was not reached, whereas the median OS was 18.10 months (95% CI = 5.10~NA) for the high-risk group ( Figure 5C). Concerning PFS, the low-risk group tended to favor favorable PFS (HR = 0.3994, 95% CI = 0.1908~0.8361, p = 0.0120, Figure 5D). For the low-risk group, the median PFS was 19.33 months (95% CI = 6.28~NA), whereas the median PFS was 2.72 months (95% CI = 1.91~24.80) for the high-risk group ( Figure 5D). Overall, the 15 lncRNAs signature appears to independently estimate OS and PFS in melanoma patients treated with anti-PD-1 monotherapy well.
Time-dependent ROC curve analysis was performed to compare the sensitivity and specificity of the prediction of OS and PFS with the 15 lncRNAs signature in the training and validation cohorts. AUC values at 12, 18, and 24 months obtained from ROC curve analysis were used to assess the prognostic accuracy. In

Response Rates Based on the lncRNAs Signature
Among advanced melanoma patients treated with anti-PD-1 monotherapy, there was a strong association between best overall response and the identified 15 lncRNAs signature. In the training cohort, 2 (7%) of 30 patients in the high-risk group versus 8 (38%) of 21 patients in the low-risk group had an objective response (CR or PR); CR were achieved by 0 (0%) of 30 patients in the high-risk group versus 3 (14%) of 21 patients in the low-risk group; 11 (37%) of 30 patients in the high-risk group had disease control (CR, PR, or SD) as their best overall response, compared to 15 (71%) of 21 patients in the low-risk group (p = 0.0236, Figure 7). In the validation cohort, 6 (32%) of 19 patients in the high-risk group versus 13 (59%) of 22 patients in the low-risk group had an objective response; CR were achieved by 1 (5%) of 19 patients in the high-risk group versus 3 (14%) of 22 patients in the low-risk group; 8 (42%) of 19 patients in the high-risk group had best overall response compared to 17 (77%) of 22 patients in the low-risk group (p = 0.0476, Figure 7).

Immune Cell Enrichment Analysis
A compendium of 39 cell types, comprising multiple adaptive and innate immune cells derived from thousands of expression profiles, was identified in the training cohort with xCell, a novel method that integrates the advantages of gene set enrichment with deconvolution approaches [35]. We identified CD4+ Th1 cells and CD8+ naïve T cells to be significantly altered (p = 0.0103 and 0.0436, respectively; Table S4).

Discussion
Over the past decades, with a deeper understanding of the pathophysiology and the manifold roles of the immune system in cancer management [39], ICI have made their way into clinics as single treatment or in multimodal settings for several tumor entities [40][41][42]. While these advances in the treatment of metastatic melanoma have improved responses and survival [43], still the majority of patients do not respond properly or respond with side effects to ICI. Although the regimens of ICI-based immunotherapy have been continuously adjusted and optimized [44][45][46], patients with melanoma still have heterogeneous ICI response, especially to anti-PD-1 monotherapy. This requires screening for the most sensitive subgroups. An early assessment, especially at the pre-treatment stage, for anti-PD-1 monotherapy responses by predictive signature is crucial for the selection of patients who are most likely to benefit from anti-PD-1 monotherapy. Hence, in this study, we combined WGCNA and LASSO to discover 15 lncRNAs predictor for patients with advanced melanoma, which can also reflect the anti-PD-1 monotherapy response, and then initially explored the potential mechanisms.
The survival curves indicate that the identified 15 lncRNAs signature distinguishes patients who are most likely to have survival benefit of anti-PD-1 monotherapy, regardless of PFS and OS. However, the signature was stronger for prediction of PFS when being validated. This is of high value, since OS might be generally influenced by non-disease specific events.
Furthermore, we used the DEGs between two consensus clusters as a starting point to get first hints why there is a survival benefit between the two clusters distinguished by the 15 lncRNAs signature by using MSigDB, which is has been seldomly performed in previous studies [47,48]. We found several immune-related cells, processes, and pathways being affected. Recent studies showed that IL2-STAT5 signaling pathway is closely related to immunity [49,50]. The activation of the MAPK8 in melanoma could trigger the massive functional natural killer cells infiltration [51]. Moreover, the enriched C2 contents jaeger metastasis dn is related to the molecular mechanisms of malignant melanoma progression and metastasis [52], and onder cdh1 targets 2 dn contributes to metastatic dissemination [53]. In C7, we also identified several immune-related cells, such as B cell, CD8 + T cell, and CD4 + T cell. Additionally, the identified pathways were also enriched for PD-L1/PD-1 axis events, which is exactly the target of anti-PD-1 therapy. In terms of adaptive immune response, CD4 + T cells are regarded as important factors regulating immune balance [54]. In the xcell analysis, we also highlighted CD4 + T cell and CD8 + T cell, which is consistent with functional analysis. T cells are prominent TILs in melanoma [55]. CD4 + T cells are associated with anti-tumor responses in melanoma [56][57][58]; CD8 + T cells also play a role herein [59]. Further, different CD8 + T cell subpopulations have predictive value in melanoma [60].
Moreover, similar to previous studies [61][62][63], we also succeeded with the identified 15 lncRNAs signature to stratify patients in the high-and low-risk groups. The high-risk group had a lower best overall response compared to the low-risk group, both in the training and in the validation cohort (Figure 7), indicating the robustness of the 15 lncRNAs signature. Therefore, one can hypothesize that the identified 15 lncRNAs signature might affect the survival benefit of advanced melanoma patients treated with anti-PD-1 monotherapy through generally immune-related cells, processes, and pathways.
The lncRNAs we identified may provide new ideas and insights for predicting the survival benefit of melanoma patients who receive anti-PD-1 monotherapy. However, one has to consider that our analyses are based on the publicly available dataset. Thus, we could not obtain all the clinic-pathological characteristics for each patient. Furthermore, we only included data of pre-treatment tumor biopsies from melanoma patients treated with anti-PD-1 monotherapy.. Although we used two completely independent data as the training cohort and the validation cohort, anti-PD-1 monotherapy in the two data is not exactly consistent. Moreover, since the population we included in this retrospective study was melanoma patients receiving anti-PD-1 monotherapy with complete follow-up data, our sample size was relatively small. Therefore, further testing and verifying of the identified 15 lncRNAs signature in prospective studies will be necessary.

Conclusions
However, we succeeded to characterize lncRNA expression profiles to identify particularly PFS benefit of melanoma patients receiving anti-PD-1 monotherapyOur analyses provide also hints that this signature affects the response of melanoma patients treated with anti-PD-1 monotherapy by influencing immune-related pathways. The 15 lncRNAs signature is therefore a novel predictor for survival in melanoma patients treated with anti-PD-1 monotherapy.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cells10050977/s1, Figure S1: The median rank and Z summary statistics of each module preservation, Figure S2 Table S1: Details of lncRNAs in the skyblue module, Table S2: Details of the DEGs, Table S3: Details of the DEGs function annotation analysis, Table S4: Immune cell enrichment analysis by xCell.