Unsupervised Analysis Reveals the Involvement of Key Immune Response Genes and the Matrisome in Resistance to BRAF and MEK Inhibitors in Melanoma

Simple Summary Drug resistance is still an imminent issue for melanoma patients even after the use of BRAFi/MEKi combination therapy in clinic. Tumor heterogeneity and heterogeneity in treatment responses make it difficult to find consensus genes and pathways in resistance to therapy. This study used an objective method to analyze published gene expression data from pre-treated tumors and drug-resistant tumors, and identified possible targets and markers for resistant tumors, which centered on the PLXNC1 gene, which promotes a pro-inflammatory tumor microenvironment. Abstract Melanoma tumors exhibit a wide range of heterogeneity in genomics even with shared mutations in the MAPK pathway, including BRAF mutations. Consistently, adaptive drug resistance to BRAF inhibitors and/or BRAF plus MEK inhibitors also exhibits a wide range of heterogeneous responses, which poses an obstacle for discovering common genes and pathways that can be used in clinic for overcoming drug resistance. This study objectively analyzed two sets of previously published tumor genomics data comparing pre-treated melanoma tumors and BRAFi- and/or MEKi-resistant tumors. Heterogeneity in response to BRAFi and BRAFi/MEKi was evident because the pre-treated tumors and resistant tumors did not exhibit a tendency of clustering together. Differentially expressed gene (DEG) analysis revealed eight genes and two related enriched signature gene sets (matrisome and matrisome-associated signature gene sets) shared by both sets of data. The matrisome was closely related to the tumor microenvironment and immune response, and five out of the eight shared genes were also related to immune response. The PLXNC1 gene links the shared gene set and the enriched signature gene sets as it presented in all analysis results. As the PLXNC1 gene was up-regulated in the resistant tumors, we validated the up-regulation of this gene in a laboratory using vemurafenib-resistant cell lines. Given its role in promoting inflammation, this study suggests that resistant tumors exhibit an inflammatory tumor microenvironment. The involvement of the matrisome and the specific set of immune genes identified in this study may provide new opportunities for developing future therapeutic methods.


Introduction
Since the discovery of BRAF as a major driver mutation in melanoma [1], inhibitors of BRAF (BRAFis) have been developed and applied effectively in clinic settings [2][3][4].Early results indicate that although BRAFis are generally effective and moderately extend life span in patients, drug resistance develops quickly.Hence, combination treatment of a BRAFi with a MEK inhibitor (MEKi) became a first-line choice of treatment, along with rapidly developed immune checkpoint inhibitors [5].
Data processing: The downloaded data was adjusted such to ensure all expression values were positive, and then the values were log2 transformed.For limma analysis, genes were ranked by a formula: "sign(logFC)*(−log10(p.Value))".Significant genes were retained by a nominal p value of <0.05 and a 1.5-fold change in either direction.A Venn diagram was drawn by Gene Venn program (https://www.bioinformatics.org/gvenn/,accessed on 10 February 2024).
Vemurafenib-resistant cells, RNA isolation, and qRT-PCR: The establishment of vemurafenib-resistant cell lines was described before [22].Briefly, the parental melanoma cell lines SK-Mel28 and 1205Lu were cultured in increasing concentrations of vemurafenib for 4-8 weeks.In the end, the cells were able to proliferate normally in a vemurafenib concentration equal to or greater than the IC50s of the parental cell lines.RNA was isolated using the PureLink™ RNA Mini Kit (Invitrogen, Carlsbad, CA, USA, cat #12183018A), and then treated with DNase before being subjected to one-step qRT-PCR.The GAPDH gene was used as an internal control (primers: 5 ′ -CCA CTCCTC CAC CTT TGA CGC -3 ′ ; 5 ′ -GAC TGAGTG TGG CAG GGACTC-3 ′ ).The primers for PLXNC1 were designed by NCBI Primer3 (forward primer: 5 ′ -GGT GGC AAT TCA TTC TGT GCT T-3 ′ , reverse primer: 5 ′ -TGA GGA TGA CAA TGG AGG CAA A-3 ′ ).Relative mRNA levels were calculated by the DDCt method, as described before [23].The experiments were repeated twice, with three replicates each time.The two-sided p values were calculated by a Student's t-test based on whether the DDCt value was different than 0 (Stata/SE18.0,College Station, TX, USA).For an initial examination of data, expression dataset GSE50509 was used for evaluation.This dataset constitutes 61 tumors from 21 patients, among which 28 were from pre-BRAFi treatment while 33 were progressed drug-resistant tumors (Table 1).Some patients had more than one tumor before treatment and/or after treatment (Table S1).First, the entire dataset was used for analysis, with an assumption that gene expression levels below the median of the normalized beadchip readings were not expressed and hence removed.Unexpectedly, the heatmap generated by the pheatmap package showed no clusters based on any of the categorical variables (Figure 1a), including pre-treated tumors vs. resistant tumors, age, sex, and response to treatment.To validate this result, PCA analysis was used to show the relative relationship of tumors and patients (Figure 1b).Similar to the heatmap results, the pre-treated tumors were mixed with progressed tumors, exhibiting no tendency of grouping by treatment.For example, one of the two progressed tumors in Patient 10 (treated with dabrafenib, Table S1) is quite close to the original two tumors from the same patient, while the other one showed a distant location.Patient 11 (treated with dabrafenib) has six tumors; two progressed tumors are closer to two of the pre-treated ones, but the other two progressed ones are distant.

3
In order to validate this result, a second dataset GSE61992 was used for a similar analysis.This dataset contained the gene expression profiles of 26 tumors from 11 patients who were treated with combined dabrafenib and trametinib (Tables 1 and S1).Again, tumors were not grouped by treatment, but they seemed to be grouped by patient in this dataset, regardless of being treated or not (Figure 1c,d).The tumors were not grouped by treatment or type of BRAF mutation.expression in GSE50509 and GSE61992, respectively.In the "response" variable, NoRes, no response; PR, partial response; RES, response; SD, stable disease.Tumors are colored by "pretreatment" and "progression" (resistant tumors).(c,d) Principal component analysis of the two datasets based on patients and treatment.EDT: early during treatment when the tumors were excised.Some patients have multiple tumors, either from pre-treatment or progression.

Analysis Using the Melanoma-Specific Gene Set Also Returned No Significant Clustering
Recognizing that the presumption above, which assumes that 50% of genes are not expressed in the dataset, is overly arbitrary, the TCGA-SKCM data was then utilized to identify all genes specifically expressed in melanoma tumors.All genes with an RSEM (RNA-Seq by Expectation-Maximization) expression level greater than one across all tumors in the TCGA-SKCM dataset were considered to be expressed in melanoma tumors.A total of 17,895 genes remained, which were termed melanoma-specific genes.Next the melanoma-specific gene sets from the GSE50509 and GSE61992 datasets were extracted, and the heatmap and PC analysis were used for clustering again.The results were identical to the first analysis, with no apparent clustering of pre-treated tumors and resistant tumors.) Heatmap of gene expression in GSE50509 and GSE61992, respectively.In the "response" variable, NoRes, no response; PR, partial response; RES, response; SD, stable disease.Tumors are colored by "pretreatment" and "progression" (resistant tumors).(c,d) Principal component analysis of the two datasets based on patients and treatment.EDT: early during treatment when the tumors were excised.Some patients have multiple tumors, either from pre-treatment or progression.

Analysis Using the Melanoma-Specific Gene Set also Returned No Significant Clustering
Recognizing that the presumption above, which assumes that 50% of genes are not expressed in the dataset, is overly arbitrary, the TCGA-SKCM data was then utilized to identify all genes specifically expressed in melanoma tumors.All genes with an RSEM (RNA-Seq by Expectation-Maximization) expression level greater than one across all tumors in the TCGA-SKCM dataset were considered to be expressed in melanoma tumors.A total of 17,895 genes remained, which were termed melanoma-specific genes.Next the melanoma-specific gene sets from the GSE50509 and GSE61992 datasets were extracted, and the heatmap and PC analysis were used for clustering again.The results were identical to the first analysis, with no apparent clustering of pre-treated tumors and resistant tumors.

Unsupervised Differential Gene Expression Analysis
Because of the above heterogeneity of gene expression profiles, differential gene expressions (DEGs) between pre-treated tumors and progressed tumors were obtained by nominal p values with a cutoff of 1.5-fold changes.Multiple comparison adjustment was not used as this would return no significant DEGs, either using 50% of all genes or using the melanoma-specific genes.
Limma analysis resulted in 380 significant genes (unadjusted p < 0.05) in the DEG set for GSE50509 and 569 significant genes (unadjusted p < 0.05) for GSE61992 (Supplementary Materials, Tables S2 and S3).A Venn diagram was used to examine the shared genes in these two sets (Figure S1).There were only 12 genes shared by these two datasets (Table 2)

Gene Set Enrichment Analysis
Each DEG set was ranked and then used for gene set enrichment analysis (GSEA).Each molecular signature set collected in MSigDB was separately used for GSEA analysis, including Hallmark gene signatures and C1 to C8 gene signatures [24][25][26].Results for GSE50509 and GSE61992 are listed in Table 3.With the adjusted p-value of 0.05 as a cutoff point, there are 26 and 22 significantly enriched gene sets for the GSE50509 and GSE61992 datasets, respectively.Interestingly only two related gene sets are found in both datasets: NABA_Matrisome and NABA_Matrisome_Associated. However, the enrichment is in two directions: the normalized enrichment scores in the GSE50509 were −2.38 and −1.83 for NABA_Matrisome and NABA_Matrisome_Associated, respectively, while that for GSE61992 were 2.35 and 2.12, respectively (Table 3, in bold).When the leading-edge components are examined, the two gene sets actually have different leading-edge genes except for one shared gene, PLXNC1, which shows the same regulation direction in two datasets.PLXNC1 was up-regulated in the resistant tumors by 1.73-and 1.51-fold in the GSE50509 and GSE61992 sets, respectively.Therefore, this shared matrisome enrichment is at least not contradictory to each other; rather, they are consistent.

Validation of PLXNC1 Over-Expression in BRAFi-Resistant Cells
To validate that PLXNC1 is indeed up-regulated in BRAFi-resistant cells, we performed a qRT-PCR experiment using two pairs of cell lines: the parental SK-Mel28 and 1205Lu cells, as well as their corresponding vemurafenib-resistant cell lines.These resistant cell lines were established in our laboratory before and reported elsewhere [22].The qRT-PCR results are shown in Figure 2. In comparison to the parental cell lines, SK-Mel28-R (SK-Mel28 resistant cell line) showed about a 4.1-fold ± 0.5 higher expression of PLXNC1 (p = 0.04) while 1205Lu-R (1205Lu resistant cells) showed a 6.2-fold ± 0.7 higher expression of PLXNC1 compared to 1205Lu cells (p = 0.02, Figure 2).

Validation of PLXNC1 Over-Expression in BRAFi-Resistant Cells
To validate that PLXNC1 is indeed up-regulated in BRAFi-resistant cells, we performed a qRT-PCR experiment using two pairs of cell lines: the parental SK-Mel28 and 1205Lu cells, as well as their corresponding vemurafenib-resistant cell lines.These resistant cell lines were established in our laboratory before and reported elsewhere [22].The qRT-PCR results are shown in Figure 2. In comparison to the parental cell lines, SK-Mel28-R (SK-Mel28 resistant cell line) showed about a 4.1-fold ± 0.5 higher expression of PLXNC1 (p = 0.04) while 1205Lu-R (1205Lu resistant cells) showed a 6.2-fold ± 0.7 higher expression of PLXNC1 compared to 1205Lu cells (p = 0.02, Figure 2).

Discussion
It is desirable to discover consensus gene regulation in tumors resistant to BRAFi and MEKi treatment for the sake of developing new treatment methods.For example, reactivation of the MAPK pathway was discovered following BRAFi treatment [9], which resulted in the combination treatment of a BRAFi and a MEKi.While there are many other similar attempts, the actual results varied with limited clinical advancement, perhaps due to the extreme heterogeneity of tumors and drug resistance responses.Such heterogeneity was reported before, but the main focus was on the MAPK pathway and melanoma-related oncogenic pathways and tumor suppressors, such as PI3K, PTEN, AKT, and CDKN2A [27], or amplification of BRAF and mutations of the RAS family [12,28].This study used published datasets and analyzed gene expression objectively and discovered that while there were a few immune-related genes showing alteration in resistant tumors, the matrisome seemed to be a significant player.
First, the heterogeneity of gene expression profiling was demonstrated by PCA analysis and a heatmap, in which the resistant tumors were mingled in with the early dissected tumors or untreated tumors, without a distinct clustering pattern.For this reason, nominal p values had to be used for DEG analysis.Only 12 genes were shared in the two datasets, among which 8 genes were in the same direction.
Five genes (IFI44, FLOR2, OLR1, PLXNC1, and SAMD9) among these eight genes were related to tumor-associated macrophages or tumor immune response.IFI44 is an interferon-a inducible protein which suppresses the immune response during viral

Discussion
It is desirable to discover consensus gene regulation in tumors resistant to BRAFi and MEKi treatment for the sake of developing new treatment methods.For example, reactivation of the MAPK pathway was discovered following BRAFi treatment [9], which resulted in the combination treatment of a BRAFi and a MEKi.While there are many other similar attempts, the actual results varied with limited clinical advancement, perhaps due to the extreme heterogeneity of tumors and drug resistance responses.Such heterogeneity was reported before, but the main focus was on the MAPK pathway and melanoma-related oncogenic pathways and tumor suppressors, such as PI3K, PTEN, AKT, and CDKN2A [27], or amplification of BRAF and mutations of the RAS family [12,28].This study used published datasets and analyzed gene expression objectively and discovered that while there were a few immune-related genes showing alteration in resistant tumors, the matrisome seemed to be a significant player.
First, the heterogeneity of gene expression profiling was demonstrated by PCA analysis and a heatmap, in which the resistant tumors were mingled in with the early dissected tumors or untreated tumors, without a distinct clustering pattern.For this reason, nominal p values had to be used for DEG analysis.Only 12 genes were shared in the two datasets, among which 8 genes were in the same direction.
Five genes (IFI44, FLOR2, OLR1, PLXNC1, and SAMD9) among these eight genes were related to tumor-associated macrophages or tumor immune response.IFI44 is an interferon-a inducible protein which suppresses the immune response during viral infection [29,30].This may be consistent with a previous report that BRAFi impacted the tumor microenvironment where the number of immune cells were reduced during BRAFi treatment [31].The FLOR2 gene is a marker for tumor-associated macrophages and is identified as an anti-inflammatory gene [32].FLOR2 is down-regulated in resistant tumors, suggesting an increased inflammatory tumor microenvironment, which is consistent with the effect of increased IFI44 gene expression in resistant tumors.As for the OLR1 gene, a recent scRNA-seq analysis in head and neck tumors showed that OLR1 was specifically expressed on tumor-associated macrophages, and was significantly associated with worse overall survival of patients [33].Similarly, PLXNC1, a member of the Plexin family composed of transmembrane domains, was also shown to be significantly associated with M2 macrophages, and indicated poor outcome in stomach adenocarcinoma and acute myeloid leukemia [34,35].The direction of gene regulation of these four genes points to poor overall survival, which is consistent with resistance to BRAFi or BRAFi/MEKi treatment.
The ARHGAP18 gene is involved in endothelial cell regulation and tumor angiogenesis [36], while the PKD2 protein which contains six transmembrane domains, plays a key role in autophagy, which is an intracellular degradation process under stressful stimuli [37].The LEF1 transcriptional factor is multi-functional, and its involvement in BRAFi/MEKi resistance has not been well established [38], but it may be related to epithelial-mesenchymal transition.
The only two shared enrichment gene sets are matrisome and matrisome-associated genes in the two GSE datasets.The matrisome is a set of gene signatures established by a proteosome method to identify extracellular matrix components, i.e., an ECM (extracellular matrix) proteome [15].Matrisome-associated genes are proteins regulating or modulating the ECM proteome [15].Together there are ~1000 matrisome and matriosome-associated genes.Both GSE datasets identified matrisome and matrisome-associated genes in the GSEA analysis, with one single shared gene in the leading gene set, PLXNC1, which was in the same direction of regulation.As the ECM or matrisome, indeed, defines the tumor microenvironment, which is closely related to immune regulation, inflammation, tumor invasiveness, and metastasis [14,16], it is no surprise that this study found that the matrisome exhibited a significant change in drug-resistant tumors.
Since PLXNC1 appeared in the shared gene set as well as in the shared enrichment gene set, its function in cellular signaling was further examined.PLXNC1 is a receptor for the SEMA7A protein, which regulates a wide range of tumor cell functions, including proliferation, invasion, migration, and angiogenesis.PLXNC1 is induced during acute inflammation and is expressed in the nervous system.Knockout of PLXNC1 led to decreased inflammatory responses and improved survival in mice [39,40].We showed that PLXNC1 was indeed up-regulated in vemurafenib-resistant SK-Mel28 and 1205Lu melanoma cells.
PLXNC1 expression was lost in metastatic melanoma as compared to matched primary tumors [41], and, therefore, PLXNC1 in melanoma is considered to be a tumor suppressor.A further piece of supporting evidence is found in a mouse model of melanoma, where PlxnC1 expression was positively correlated with Ednrb, suggesting a role in the Ednrb/EDNRBmediated suppressive effect on melanoma development [42].However, in hepatocellular carcinoma and stomach cancer, PLXNC1 promotes cancer cell proliferation and metastasis, and is associated with overall poorer survival [43,44].In stomach cancer, PLXNC1 was significantly associated with M2 macrophages and poor outcome [34].In the context of the BRAFi and MEKi treatment condition, increased expression of PLXNC1 perhaps indicated an increased inflammatory tumor microenvironment.
It was reported that BRAFi can create a pro-inflammatory tumor microenvironment via reduced recruitment of immunosuppressive cells and increased immune effector cells, as well as an increased ability of immune effector cells to trigger cancer killing [45].This inflammatory tumor microenvironment could be subjective to dynamic regulation when drug resistance develops [45].
As the overall sample size is not large enough in this study, the resulting common genes and signature gene sets need further validation with larger datasets, or validation in a laboratory.If results are validated, it is possible to develop PLXNC1 as a progression biomarker for resistant tumors, or as a future treatment target.

Conclusions
Consistent with the functions of the eight shared genes found in the two datasets, the matrisome changes with tumor progression and is also connected to the tumor immune response [14,16], and is a possible target for developing targeted therapy.Therefore, overall, this study reveals previously unknown aspects of the BRAFi and MEKi resistance phenotype, i.e., modulation of the matrisome and a few specific immune-related genes that were not identified in previous studies.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers16132313/s1, Figure S1: Venn Diagram for shared genes in the two DEG sets.Table S1: Details about patients and tumors.Table S2: Significant genes in GSE50509.Table S3: Significant genes in GSE61992.Institutional Review Board Statement: Ethical review and approval were waived for this study due to no involvement of personal identification information.
Informed Consent Statement: Not applicable.

. 1 .
The Clustering of Tumors and Patients Reveals No Clustering by Pre-Treatment or Post-Treatment, or by Any Other Criteria

Figure 1 .
Figure 1.Melanoma tumors and patients are not grouped by drug resistance.(a,b) Heatmap of geneexpression in GSE50509 and GSE61992, respectively.In the "response" variable, NoRes, no response; PR, partial response; RES, response; SD, stable disease.Tumors are colored by "pretreatment" and "progression" (resistant tumors).(c,d) Principal component analysis of the two datasets based on patients and treatment.EDT: early during treatment when the tumors were excised.Some patients have multiple tumors, either from pre-treatment or progression.

Figure 1 .
Figure 1.Melanoma tumors and patients are not grouped by drug resistance.(a,b) Heatmap of gene expression in GSE50509 and GSE61992, respectively.In the "response" variable, NoRes, no response; PR, partial response; RES, response; SD, stable disease.Tumors are colored by "pretreatment" and "progression" (resistant tumors).(c,d)Principal component analysis of the two datasets based on patients and treatment.EDT: early during treatment when the tumors were excised.Some patients have multiple tumors, either from pre-treatment or progression.

Figure 2 .
Figure 2. PLXNC1 is up-regulated at an mRNA level in BRAFi-resistant melanoma cell lines.The qRT-PCR experiment was performed with GAPDH as an internal control, and the DDCt method was used to calculate the relative expression levels of the resistant cells normalized to the parental control cells.

Figure 2 .
Figure 2. PLXNC1 is up-regulated at an mRNA level in BRAFi-resistant melanoma cell lines.The qRT-PCR experiment was performed with GAPDH as an internal control, and the DDCt method was used to calculate the relative expression levels of the resistant cells normalized to the parental control cells.

Table 1 .
Tumor characteristics from the two datasets.

Table 2 .
The 12 shared genes in the two datasets.

Table 3 .
The GSEA results.Bolded: shared enriched pathways by the two datasets.