An estrogen response-related signature predicts response to immunotherapy in melanoma

Background: Estrogen/estrogen receptor signaling influences the tumor microenvironment and affects the efficacy of immunotherapy in some tumors, including melanoma. This study aimed to construct an estrogen response-related gene signature for predicting response to immunotherapy in melanoma.

Methods: RNA sequencing data of 4 immunotherapy-treated melanoma datasets and TCGA melanoma was obtained from open access repository. Differential expression analysis and pathway analysis were performed between immunotherapy responders and non-responders. Using dataset GSE91061 as the training group, a multivariate logistic regression model was built from estrogen response-related differential expression genes to predict the response to immunotherapy. The other 3 datasets of immunotherapy-treated melanoma were used as the validation group. The correlation was also examined between the prediction score from the model and immune cell infiltration estimated by xCell in the immunotherapy-treated and TCGA melanoma cases.

Results: “Hallmark Estrogen Response Late” was significantly downregulated in immunotherapy responders. 11 estrogen response-related genes were significantly differentially expressed between immunotherapy responders and non-responders, and were included in the multivariate logistic regression model. The AUC was 0.888 in the training group and 0.654–0.720 in the validation group. A higher 11-gene signature score was significantly correlated to increased infiltration of CD8+ T cells (rho=0.32, p=0.02). TCGA melanoma with a high signature score showed a significantly higher proportion of immune-enriched/fibrotic and immune-enriched/non-fibrotic microenvironment subtypes (p<0.001)–subtypes with better response to immunotherapy–and significantly better progression-free interval (p=0.021).

Conclusion: In this study, we identified and verified an 11-gene signature that could predict response to immunotherapy in melanoma and was correlated with tumor-infiltrating lymphocytes. Our study suggests targeting estrogen-related pathways may serve as a combination strategy for immunotherapy in melanoma.

Introduction

Melanoma is an aggressive malignant skin cancer that causes the majority of deaths from skin cancers. Chemotherapies by multiple target therapies have been developed during the last decades to combat molecular defects of melanoma, including BRAF inhibitors vemurafenib and dabrafenib (1). However, although these drugs are highly effective for patients with ^V600BRAF-mutated melanomas, which account for approximately half of metastatic melanomas, many patients develop resistance within a relatively short period (2). Melanomas are among the most immunogenic tumors therefore studies of immunotherapy (mostly immune checkpoint blockade [ICB] therapy) for metastatic melanoma have received considerable attention. Immune checkpoint inhibitors such as anti-programmed cell death 1 (anti-PD-1) significantly improve relapse-free survival in patients with resected stage III/IV melanomas (3, 4). However, less than 50% of patients could get tumor regression and long-term durable cancer control under ICB therapy, and chronic immune-related adverse events appear to be more common after ICB therapy (5, 6). Therefore, it is important to identify patients who can benefit from ICB therapy and find new strategies to enhance its effectiveness.

Gender influences the progression of melanoma during all phases, with women showing a lower incidence and lower risk of lymph node invasion and visceral metastases compared to men (7). Clinical data showed that female patients with advanced melanoma may not benefit as much from combination ICB treatment as male patients and estrogen level may serve as an important biomarker associated with ICB therapy response (8). Estrogen/estrogen receptor (ER) signaling influences the tumor microenvironment (TME) and affects the efficacy of ICB treatment in certain tumors (9, 10). The immune cells in the melanoma TME have complex crosstalk with the tumor cells which affects the response to treatments (11). A recent study by Chakraborty et al. highlighted that inhibition of estrogen signaling influences intratumoral macrophage polarization in melanoma, increasing ICB efficacy (12). A former study in lung cancer cells identified the selective estrogen receptor degrader fulvestrant as the top compound that increased tumor sensitivity to immune-mediated lysis (13). All these studies suggest estrogen plays a role in antitumor immunity, and an estrogen response-related signature may predict the response to immunotherapy in melanoma.

The objective of the present study was to construct an estrogen response-related gene signature and evaluate its predictive ability for immunotherapy response in melanoma. Additionally, this study investigated the feasibility of combining endocrine therapy with immunotherapy in melanoma.

Materials and methods

Patient data acquisition

RNA sequencing (RNA-seq) data of immunotherapy-treated melanoma were obtained from Riaz et al. (anti-PD-1, N=51, https://github.com/riazn/bms038, Gene Expression Omnibus [GEO] accession number GSE91061) (14), Lauss et al. (adoptive T-cell therapy using tumor-infiltrating lymphocytes, N=25, GEO accession number GSE100797) (15), Hugo et al. (anti-PD-1, N=27, GEO accession number GSE78220) (16), and Van Allen et al. (anti-CTLA4, N=39, www.cbioportal.com, study id skcm_dfci_2015) (17). TCGA melanoma RNA expression data (transcript per million) were downloaded from GEO (accession number GSE62944) (18). Survival data of TCGA were retrieved from Liu et al. (19).

Differential expression analysis

Differential expression analysis was performed using R package DESeq2 (version 1.30.0) (20) using gene raw counts data from https://github.com/riazn/bms038. Genes with absolute fold change>1.5 and adjusted p-value (FDR) <0.05 were selected as differentially expressed genes. The list of estrogen response-related genes (N=299) was obtained by combining the genes from the Molecular Signatures Database (MSigDB) “Hallmark Estrogen Response Early” and “Hallmark Estrogen Response Late”.

Pathway analysis

50 Hallmark gene sets were downloaded from the Molecular Signature Database (MsigDB, version 7.5.1). Gene set enrichment analysis (GSEA) was performed using R package clusterProfiler (version 3.18.0) (21). Default settings in GSEA function were used excepted the following parameter: eps=0, seed=12345, and pvalueCutoff=1. Pathways with adjusted p-value <0.05 and normalized enrichment score (NES)>1 or NES<-1 were defined as significantly upregulated or downregulated pathways, respectively.

Immune cell infiltration analysis

XCell and MCPcounter immune cell enrichment scores were estimated using R package xCell (version 1.1.0) (22) and MCPcounter (version 1.2.0) (23), and gene expression data in transcript per million (TPM) were used as input. Fragments per kilobase per million mapped fragments (FPKM) data of GSE91061 was downloaded from GEO (accession number GSE91061) and transformed to TPM in R. TCGA microenvironment subtypes (immune-enriched/fibrotic [IE/F], immune-enriched/non-fibrotic [IE], fibrotic [F], and desert [D]) were downloaded from Bagaev et al. (n=463) (24). H&E image-based tumor-infiltrating lymphocyte (TIL) percentage estimation was obtained from Saltz et al. (n=383) (25).

Statistical analysis

All statistical analyses were conducted using R software (version 4.0.3). R package stats (version 4.0.3) was used for univariate and multivariate logistic regression. Gene expression data (log2TPM) of the 4 immunotherapy-treated melanoma datasets and TCGA melanoma cases were first normalized by z-score normalization separately. The immunotherapy response data and the expression levels (z-score) of the 11 estrogen response-related genes from Riaz et al. (14) (GSE91061) were used as the input for the construction of the 11-gene prediction model. The model was then applied to data from Lauss et al. (15), Hugo et al. (16), Van Allen et al. (17) and TCGA (18) to calculate the prediction score. R package “pROC”(version 1.18.0) was used to plot the receiver operating characteristic (ROC) curve and calibration curve (26). Wilcoxon rank-sum test and Chi-squared test were used in the comparison between two groups for continuous and categorical variables, respectively. Log-rank test was used for two-group survival comparison in Kaplan-Meier plot.

Results

Estrogen response signatures were downregulated in ICB responders.

Using the pre-treatment RNA-seq data of 51 ICB-treated melanoma (GSE91061), we first performed differential expression analysis between ICB responders (N=10) and non-responders (N=41) (Supplementary Table 1). 77 upregulated genes and 155 downregulated genes (fold change>1.5 and adjusted p-value<0.05) were identified in ICB responders as compared to non-responders (Figure 1A, Supplementary Table 2, Supplementary Figure 1). To identify pathways changed between ICB responders and non-responders, we applied gene set enrichment analysis (GSEA) using the 50 hallmark gene sets. Immune-related pathways such as “Hallmark Allograft Rejection”, “Hallmark Interferon Gamma Response”, “Hallmark IL6 JAK STAT3 Signaling”, “Hallmark Inflammatory Response”, and “Hallmark IL2 STAT5 Signaling” were significantly activated in ICB responders (Figure 1B, Supplementary Table 3). 5 pathways were significantly downregulated in ICB responders (Figure 1C, Supplementary Table 3), among which “Hallmark Estrogen Response Late” was an estrogen response-related pathway.

FIGURE 1

Figure 1 Estrogen response signatures were activated in ICB responders. (A) Differentially expressed genes between ICB responders (n=10) and non-responders (n=41). Upregulated genes and downregulated genes (fold change>1.5 and adjusted p-value<0.05) were marked in red and blue, respectively. 11 estrogen response-related genes were labeled. (B, C) Top 10 upregulated (B) and downregulated (C) pathways in ICB responders as compared to non-responders. 50 hallmark gene sets were analyzed using GSEA. Gene sets with adjusted p-value<0.05 were labeled in grey. NES, normalized enrichment score.

Establishment and evaluation of an 11-gene estrogen response-related signature to predict immunotherapy response

Among the 232 differentially expressed genes between ICB responders and non-responders, 11 genes (AGR2, KLK11, PKP3, ELF3, FGFR3, TRIM29, SFN, KLK10, SCNN1A, CA12, and ESRP2, Figure 1A) were estrogen response-related genes. Univariate and multivariate logistic regression showed that none of the 11 genes was significantly associated with ICB response in the 51 patients (Supplementary Table 4). Considering estrogen response is a complex pathway involving many genes, we included all 11 genes and developed a multivariate logistic regression model for ICB response in melanoma. The AUC of the prediction model was 0.888 (95% confidence interval 0.786-0.990, by 500 times bootstrap resampling, Figure 2A). The calibration curve showed relatively high agreement between the predicted and the actual observation of the ICB responder in the patients with a high predicted probability of being an ICB responder (Figure 2B). Besides, heatmap clustering showed the 11-gene signature score could reflect the differentially expressed gene patterns between ICB responders and non-responders (Supplementary Figure 2). We further tested the prediction ability of this model using pre-treated RNA expression data of immunotherapy-treated melanoma from Lauss et al. (adoptive T-cell therapy, n=25) (15), Hugo et al. (anti-PD-1, n=27) (16), and Van Allen et al. (anti-CTLA4, n=39) (17) of which the AUCs were 0.720, 0.654 and 0.692 (Figures 2C–E), respectively.

FIGURE 2

Figure 2 Evaluation of the 11-gene estrogen response-related signature. (A) Receiver operating characteristic (ROC) curve of the prediction model in GSE91061 (n=51). (B) Calibration curve for the logistic regression model in GSE91061. Dashed line, prediction calibration curve. Solid line, standard curve. (C-E) ROC curves of the prediction model in Lauss et al. (n=25) (C) (15), Hugo et al. (n=27) (D) (16), and Van Allen et al. (n=39) (E) (17).

High 11-gene signature score was related to high level of CD8+ T cells and good prognosis

Estrogen response was reported to play an important role in the tumor microenvironment (TME) (27). To study the pathways related to the 11-gene signature score, we performed DE analysis using the median 11-gene signature score as the cutoff for the 51 patients in GSE91061. GSEA showed that “Hallmark Interferon Gamma Response”, “Hallmark Allograft Rejection”, “Hallmark Myogenesis”, “Hallmark Interferon Alpha Response”, and “Hallmark Inflammatory Response” were significantly upregulated in the high-signature score group (Supplementary Table 5). As expected, “Hallmark Estrogen Response Late” was downregulated in the high-signature score group (Supplementary Table 5). We further tested the correlations of the 11-gene signature score to xCell and MCPcounter immune cell enrichment scores. We found that the 11-gene signature score was significantly correlated with xCell scores of B cells (rho=0.32, p=0.021), CD4 memory T cells (rho=0.33, p=0.017), and CD8+ T cells (rho=0.32, p=0.02) (Figure 3A, Supplementary Table 6), while no significant correlation was identified between the 11-gene signature score and MCPcounter scores (Supplementary Figure 3A). As GSE91061 dataset has a limited number of patients, we then evaluated the relations of the 11-gene signature score to immune cells in TCGA melanoma cases (n=469). Patients with high 11-gene signature scores had significantly increased xCell scores of B cells (p<0.001), CD4+ memory T-cells (p<0.001), M1 macrophages (p=0.004), macrophages (p=0.001), and Tregs (p=0.0076), and a trend of high CD8+ T-cells score (p=0.091) (Figure 3B, Supplementary Table 7). Using MCPcounter scores, we also found higher scores of B lineage (p<0.05), monocytic lineage (p<0.05), CD8+ T cells (p=0.07), and cytotoxic lymphocytes (p=0.06) in TCGA high 11-gene signature score group (Supplementary Figure 3B, Supplementary Table 7). A Higher TIL percentage (p=0.078) was also identified in TCGA high 11-gene signature score group using H&E image-based TIL data (Supplementary Figure 3C) (25). Bagaev et al. (24) identified four microenvironment subtypes (immune-enriched/fibrotic, immune-enriched/non-fibrotic, fibrotic, and desert) in TCGA patients, among which the immune-enriched/fibrotic and immune-enriched/non-fibrotic subtypes had a higher level of immune activation and better response to ICB as compared to the other two subtypes. We also observed a significantly higher proportion of immune-enriched/fibrotic and immune-enriched/non-fibrotic subtypes (49.3% vs 33.8%, chi-squared test, p-value<0.001) in the TCGA melanoma cases with high 11-gene signature score (Figure 3C). Besides, TCGA melanoma cases with a high 11-gene signature score had significantly better progression-free interval (p=0.021) and a tendency towards better overall survival (p=0.055) (Figure 3D).

FIGURE 3

Figure 3 High 11-gene signature score was related to high level of CD8+ T cells and good prognosis. (A) Correlations between the 11-gene signature score and the xCell immune cell enrichment scores in GSE91061 (n=51). Xcell immune cell enrichment scores for different types of immune cells were evaluated using R package xCell. Two-sided spearman’s correlation test. A liner regression line was added to help better visualize the correlation using R package ggplot2 for plots with spearman’s correlation test p-values less than 0.05. (B) The xCell immune cell enrichment scores for different types of immune cells in TCGA melanoma cases (n=469) with high (n=234) and low (n=235) 11-gene signature scores. Median was used as the cutoff for groups with high and low signature scores. Two-sided Wilcoxon rank-sum test. (C) Microenvironment subtypes of TCGA melanoma cases with high (n=229) and low (n=234) 11-gene signature scores (median as cutoff). Four microenvironment subtypes (immune-enriched/fibrotic [IE/F], immune-enriched/non-fibrotic [IE], fibrotic [F], and desert [D]) were downloaded from Bagaev et al. (24). 463 of the 469 melanoma cases in (B) had microenvironment subtypes. Percentages of each subtype in the high and low signature score groups were labeled. Chi-squared test. (D) Overall survival and progression-free interval of TCGA melanoma cases with high and low 11-gene signature scores (median as cutoff). Log-rank test.

Inhibition of estrogen signaling may contribute to better ICB response in melanoma

To identify potential drugs that could turn melanoma from ICB non-responders to ICB responders, we queried Connectivity map (CMap) Touchstone datasets for drugs that induced similar gene expression patterns as the DE gene expression profile between ICB responders and ICB non-responders in GSE91061. In melanoma cell line-A375, we found that an estrogen receptor antagonist, Y-134 was among the top 20 drugs which had the most similar (ranked by normalized connectivity score, FDR<0.05) induced gene expression pattern as the DE genes described above (Figure 4A, Supplementary Table 8), and an estrogen agonist, DY-131 had negative connectivity score (FDR<0.05), which indicates the induced gene expression pattern of this drug was opposing to our input DE gene expression profile (Supplementary Table 8). Further GSEA using previously published gene signature by selective estrogen receptor modulator (SERM) or selective estrogen receptor degrader (SERD) (28) showed that downregulated gene signature by SERM or SERD was enriched in melanoma ICB-responders (Figure 4B, Supplementary Figure 4, Supplementary Table 9). All these data suggesting combining ICB with the inhibition of estrogen signaling may lead to improved ICB response in melanoma.

FIGURE 4

Figure 4 Inhibition of estrogen signaling may contribute to better ICB response in melanoma. (A) Top 20 candidate drugs which may turn melanoma from ICB non-responder to responder. Connectivity map (CMap) Touchstone datasets was queried to identify drugs that could induce similar gene expression pattern as the DE gene expression profile between ICB responders and ICB non-responders in GSE91061. Adjusted p-values<0.05 for all drugs. (B) Downregulated gene signature by SERM or SERD was enriched in melanoma ICI-responders. GSEA analysis was performed in GSE91061 using signature “Frasor Response to SERM or Fulvestrant DN (downregulated genes)”.

Discussion

Melanoma accounts for the majority of deaths from skin cancer. For non-resectable/metastatic melanoma, much progress has been made in targeted therapy and immunotherapy. However, there is still a lack of practical prognostic markers. Estrogen receptors are broadly expressed in many cell types involved in the innate and adaptive immune responses and play a potential immune regulatory role in the TME (29). A recent study elucidated estrogen signaling influences intratumoral macrophage polarization in melanoma and consequently enhances ICB therapy efficacy (12). These findings suggest that, estrogen response-related genes can potentially be predictive markers for ICB therapy response in melanoma. In the current study, using estrogen response-related genes that were differentially expressed between ICB responders and non-responders in melanoma, we constructed an 11-gene immunotherapy response prediction signature with stable predictive performance in different melanoma datasets treated with multiple types of immunotherapies (anti-PD-1, anti-CTLA-4 and adoptive T-cell therapy). This signature was also significantly correlated with the infiltration of multiple types of tumor immune cells and was prognostic for overall survival in melanoma.

Estrogen and estrogen response-related genes play a modulatory role in melanoma progression, probably through influencing antitumor immunity. Melanoma incidence has a gender divergence. Slightly higher rates of melanoma have been reported for young women (20-45 years) which subsequently decrease after 45 years of age (30). On the contrary, melanoma incidence progressively increases in males after 50 years of age (31). An increased risk of melanoma was associated with early age at menarche and late age at menopause, which highlighted the role of sex steroid hormones in melanoma (32). Estrogen signaling accelerates the progression of different estrogen-insensitive tumor models by contributing to deregulated myelopoiesis by both driving the mobilization of myeloid-derived suppressor cells (MDSC) and enhancing their intrinsic immunosuppressive activity (33). A study by Chakraborty et al. revealed that estrogen signaling affects intratumoral macrophage polarization in melanoma and increased ICB efficacy (12). ERβ has been reported to be the predominant ER subtype in melanoma and could represent a marker for metastatic potential and prognosis (34). ERβ activation might impair melanoma development through the inhibition of the PI3K/Akt pathway and displays a protective role in the metastatic process (34, 35). In this study, GSEA revealed that estrogen response-related pathway “Hallmark Estrogen Response Late” was significantly downregulated in ICB responders. These findings indicate the potential of estrogen response-related gene signature in predicting ICB therapy response in melanoma.

The prognostic model proposed in the present study was composed of 11 estrogen response-related genes (AGR2, CA12, ELF3, ESRP2, FGFR3, KLK10, KLK11, PKP3, SFN, SCNN1A, TRIM29), among which many genes correlated with tumor immunity and immunotherapy. Carbonic anhydrase 12 (CA12) mediated the survival of macrophages in relatively acidic TME, while on the other hand, it induced macrophage production of large amounts of C-C motif chemokine ligand 8 (CCL8), which enhanced cancer cell epithelial-mesenchymal transition and facilitated tumor metastasis (36). CA12 was included in a former gene signature for predicting the prognosis of uveal melanoma (37). Fibroblast growth factor receptor 3 (FGFR3) is among the receptor tyrosine kinases which may be activated via autocrine circuits in melanoma (38). FGFR3 is a biomarker of immune infiltration and immunotherapy response of bladder cancer (39). Kallikrein-related peptidases (KLKs) have been reported to possess novel functions in innate immunity and inflammation. KLK10 is expressed in the follicular dendritic cells that are essential for the maturation of B cells (40). KLK10 was dynamically regulated in T cells in vitro in response to viral antigens and in activated monocytes, pointing to its activities in the development of adaptive and innate immune function (41). Plakophilin 3 (PKP3) encodes a component of desmosomes with mechanical barrier function in the skin and other normal tissues. PKP3 expression was revealed to be markedly elevated in melanoma metastasis lacking immune gene signature and was strongly associated with decreased patient survival (42). SCNN1A encodes the α subunit of epithelial sodium channel and was reported to be relevant to tumor progression in a variety of cancers (43). Lou et al. revealed that SCNN1A involves in tumor immune process by influencing tumor immune cell infiltration (44). TRIM29 is a negative regulator of NK cell functions (45). TRIM29 expression was higher in patients with higher TIL and proved to be related to immune dysfunction in colorectal cancer (46). In this study, the 11-gene signature score showed a correlation with tumor immune infiltration, which may be the basis of its predictiveness in the response to different immunotherapies.

In this study, we found that estrogen receptor antagonist Y-134 could induce a similar gene expression pattern as the DE gene expression profile between ICB responders and non-responders while estrogen agonist DY-131 could induce an opposite gene expression pattern. This discovery suggests the possibility of using anti-estrogen therapy to enhance the efficacy of immunotherapy. Accumulating evidence from experimental and clinical studies has revealed the multifaceted immunomodulatory effects of endocrine therapies, especially in the modulation of TME. SERM like tamoxifen and raloxifene could affect the functional differentiation and immunostimulatory capacity of dendritic cells (47). CARMINA 02 trial assessed 86 pre- and post-neoadjuvant endocrine therapy (NET) tumor samples and found significantly increased TIL numbers in post-NET samples of responders (48). Preclinical studies indicated that SERD, a class of ERα antagonists, interacts with ER-positive immune cells in the TME such as MDSCs, TILs, and other selected immune cell subpopulations. SERD-induced inhibition of MDSCs and concurrent actions on CD8+ and CD4+ T-cells promote the interaction of immune checkpoint inhibitors with breast cancer cells and augment the curative effect (10). In melanoma, inhibition of estrogen signaling affects intratumoral macrophage polarization and increased ICB efficacy (12). Our GSEA using previously published gene signatures by SERM or SERD also showed that downregulated gene signature by SERM or SERD was enriched in melanoma ICB-responders. These findings raise the possibility of using anti-estrogens as an approach to enhance the effectiveness of ICB therapies in melanoma, but further research is needed.

There are several limitations to this study. Firstly, the datasets in this study were all retrospective data of small sample sizes from public databases. More prospective real-world data are needed to verify the clinical utility of this model. Secondly, the conclusions obtained were based on bioinformatics analysis and require further validation in vivo and in vitro. Thirdly, the use of solely estrogen response-related genes to build a prognostic model may have excluded other prognostic genes in melanoma, which limits the overall predictive capacity of the model. Finally, this model did not incorporate other risk factors and clinical parameters that may impact prognosis, which may lead to underestimating the true prognostic capacity of the model.

In conclusion, in this study, we identified and verified an 11-gene signature that could predict response to immunotherapy in melanoma. This model was a prognostic factor for melanoma patients and was correlated with TIL. The results of this study provide new perspectives regarding the potential strategies for combining immunotherapy with endocrine therapy for treating melanoma patients.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE91061 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE100797 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78220 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944 http://www.cbioportal.org/study/summary?id=skcm_dfci_2015.

Author contributions

TD and ML designed the experiment. TD and ML analyzed and interpreted the data. X-FT, LC, YL and Y-FZ assisted with the data and figure integration. TD, X-FT provided methodology, reviewing, and editing. WZ and J-HZ conceived and supervised the study. ML wrote the manuscript with all authors' contributions to writing and providing feedback. WZ and J-HZ acted as guarantors and corresponding authors for this study. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 82202179).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1109300/full#supplementary-material

Supplementary Figure 1 | Schematic showing the construction and evaluation of the 11-gene estrogen response related ICB response prediction signature.

Supplementary Figure 2 | Heatmap of the differentially expressed genes showing clustering of the 11-gene signature score high and low groups. Z-scores of differentially expressed genes from were used. Heatmap was plotted using R function heatmap.3. Canberra distance between samples was calculated and hclust=ward was used in the hierarchical clustering. ICB responders and non-responders were labeled with black and grey. Patients with high and low 11-gene signature scores (median as cutoff) were labeled with red and blue. The 11 genes used in the final signature were labeled in yellow.

Supplementary Figure 3 | (A) Correlations between the 11-gene signature score and the MCPcounter immune cell enrichment scores in GSE91061 (n=51). MCPcounter immune cell enrichment scores for different types of immune cells were evaluated using R package MCPcounter. Two-sided spearman’s correlation test. A liner regression line with a confidence interval of 0.95 was added to help better visualize the correlation using R package ggplot2. Of note, the slope of the calculated line by ggplot2 “geom_smooth” function (method=“lm”) may not be perfectly matched to the spearman’s rank correlation coefficient (rho). (B) The MCPcounter immune cell enrichment scores for different types of immune cells in TCGA melanoma cases (n=469) with high (n=234) and low (n=235) 11-gene signature scores. Median was used as the cutoff for groups with high and low signature scores. Two-sided Wilcoxon rank-sum test. (C) H&E image-based tumor-infiltrating lymphocyte (TIL) percentage of TCGA melanoma cases with high and low 11-gene signature scores (median as cutoff, high, n=199; low, n=184).

Supplementary Figure 4 | GSEA plot of “Frasor Response to SERM or Fulvestrant Up (upregulated genes)”.

Supplementary Table 1 | Clinical characteristics of melanoma cases used in the study.

Supplementary Table 2 | List of differentially expressed genes between ICB responders and non-responder in GSE91061.

Supplementary Table 3 | GSEA between ICB responders and non-responders in GSE91061 using 50 hallmark gene sets.

Supplementary Table 4 | Univariate and multivariate analyses of the 11-estrogen response-related genes in GSE91061.

Supplementary Table 5 | GSEA between high and low 11-gene signature score groups in GSE91061 using 50 hallmark gene sets.

Supplementary Table 6 | Correlations between the 11-gene signature score and the xCell immune cell enrichment scores in GSE91061.

Supplementary Table 7 | Correlations between the 11-gene signature score and the xCell and MCPcounter immune cell enrichment scores in TCGA melanoma cases.

Supplementary Table 8 | List of drugs that could induce significant (adjusted p.value<0.05) similar (or inverse) gene expression pattern as the DE gene expression profile between ICB responders and ICB non-responders in GSE91061.

Supplementary Table 9 | List of genes in “Frasor Response to SERM or Fulvestrant Up” and “Frasor Response to SERM or Fulvestrant Down”.

Abbreviations

ICB, Immune checkpoint blockade; PD-1, Programmed cell death 1; ER, Estrogen receptor; TME, Tumor microenvironment; GSEA, Gene set enrichment analysis; NES, Normalized enrichment score; TPM, Transcript per million; TIL, Tumor-infiltrating lymphocyte; ROC, Receiver operating characteristic; SERM, Selective estrogen receptor modulator; SERD, Selective estrogen receptor degrader; MDSC, Myeloid-derived suppressor cells; CA12, Carbonic anhydrase 12; CCL8, C-C motif chemokine ligand; FGFR3, Fibroblast growth factor receptor 3; KLK, Kallikrein-related peptidase; PKP3, Plakophilin 3; NET, Neoadjuvant endocrine therapy.

References

1. Robert C, Grob JJ, Stroyakovskiy D, Karaszewska B, Hauschild A, Levchenko E, et al. Five-year outcomes with dabrafenib plus trametinib in metastatic melanoma. N Engl J Med (2019) 381(7):626–36. doi: 10.1056/NEJMoa1904059

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Shi H, Hugo W, Kong X, Hong A, Koya RC, Moriceau G, et al. Acquired resistance and clonal evolution in melanoma during BRAF inhibitor therapy. Cancer Discovery (2014) 4(1):80–93. doi: 10.1158/2159-8290.CD-13-0642

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Eggermont AMM, Blank CU, Mandala M, Long GV, Atkinson V, Dalle S, et al. Adjuvant pembrolizumab versus placebo in resected stage III melanoma. N Engl J Med (2018) 378(19):1789–801. doi: 10.1056/NEJMoa1802357

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Patrinely JR, Young AC, Quach H, Williams GR, Ye F, Fan R, et al. Survivorship in immune therapy: assessing toxicities, body composition and health-related quality of life among long-term survivors treated with antibodies to programmed death-1 receptor and its ligand. Eur J Cancer. (2020) 135:211–20. doi: 10.1016/j.ejca.2020.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Patrinely JR, Johnson R, Lawless AR, Bhave P, Sawyers A, Dimitrova M, et al. Chronic immune-related adverse events following adjuvant anti-PD-1 therapy for high-risk resected melanoma. JAMA Oncol (2021) 7(5):744–8. doi: 10.1001/jamaoncol.2021.0051

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Carlino MS, Larkin J, Long GV. Immune checkpoint inhibitors in melanoma. Lancet (2021) 398(10304):1002–14. doi: 10.1016/S0140-6736(21)01206-X

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Joosse A, de Vries E, Eckel R, Nijsten T, Eggermont AMM, Hölzel D, et al. Gender differences in melanoma survival: female patients have a decreased risk of metastasis. J Invest Dermatol (2011) 131(3):719–26. doi: 10.1038/jid.2010.354

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Jang SR, Nikita N, Banks J, Keith SW, Johnson JM, Wilson M, et al. Association between sex and immune checkpoint inhibitor outcomes for patients with melanoma. JAMA Netw Open (2021) 4(12):e2136823. doi: 10.1001/jamanetworkopen.2021.36823

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wang T, Jin J, Qian C, Lou J, Lin J, Xu A, et al. Estrogen/ER in anti-tumor immunity regulation to tumor cell and tumor microenvironment. Cancer Cell Int (2021) 21(1):295. doi: 10.1186/s12935-021-02003-w

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Márquez-Garbán DC, Deng G, Comin-Anduix B, Garcia AJ, Xing Y, Chen H-W, et al. Antiestrogens in combination with immune checkpoint inhibitors in breast cancer immunotherapy. J Steroid Biochem Mol Biol (2019) 193:105415. doi: 10.1016/j.jsbmb.2019.105415

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Marzagalli M, Ebelt ND, Manuel ER. Unraveling the crosstalk between melanoma and immune cells in the tumor microenvironment. Semin Cancer Biol (2019) 59:236–50. doi: 10.1016/j.semcancer.2019.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Chakraborty B, Byemerwa J, Shepherd J, Haines CN, Baldi R, Gong W, et al. Inhibition of estrogen signaling in myeloid cells increases tumor immunity in melanoma. J Clin Invest (2021) 131(23). doi: 10.1172/JCI151347

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hamilton DH, Griner LM, Keller JM, Hu X, Southall N, Marugan J, et al. Targeting estrogen receptor signaling with fulvestrant enhances immune and chemotherapy-mediated cytotoxicity of human lung cancer. Clin Cancer Res (2016) 22(24):6204–16. doi: 10.1158/1078-0432.CCR-15-3059

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Riaz N, Havel JJ, Makarov V, Desrichard A, Urba WJ, Sims JS, et al. Tumor and microenvironment evolution during immunotherapy with nivolumab. Cell (2017) 171(4):934–49.e16. doi: 10.1016/j.cell.2017.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat Commun (2017) 8(1):1–11. doi: 10.1038/s41467-017-01460-0

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell (2016) 165(1):35–44. doi: 10.1016/j.cell.2016.02.065

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science (2015) 350(6257):207–11. doi: 10.1126/science.aad0095

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Rahman M, Jackson LK, Johnson WE, Li DY, Bild AH, Piccolo SR. Alternative preprocessing of RNA-sequencing data in the cancer genome atlas leads to improved analysis results. Bioinformatics (2015) 31(22):3666–72. doi: 10.1093/bioinformatics/btv377

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell (2018) 173(2):400–16.e11. doi: 10.1016/j.cell.2018.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18(1):1–14. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol (2016) 17(1):1–20. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell (2021) 39(6):845–65.e7. doi: 10.1016/j.ccell.2021.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Saltz J, Gupta R, Hou L, Kurc T, Singh P, Nguyen V, et al. Spatial organization and molecular correlation of tumor-infiltrating lymphocytes using deep learning on pathology images. Cell Rep (2018) 23(1):181–93.e7. doi: 10.1016/j.celrep.2018.03.086

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez J-C, et al. pROC: an open-source package for r and s+ to analyze and compare ROC curves. BMC Bioinf (2011) 12(1):1–8. doi: 10.1186/1471-2105-12-77

CrossRef Full Text | Google Scholar

27. Chakraborty B, Byemerwa J, Krebs T, Lim F, Chang C-Y, McDonnell DP. Estrogen receptor signaling in the immune system. Endocrine Rev (2023) 44(1):117–41. doi: 10.1210/endrev/bnac017

CrossRef Full Text | Google Scholar

28. Frasor J, Danes JM, Komm B, Chang KC, Lyttle CR, Katzenellenbogen BS. Profiling of estrogen up-and down-regulated gene expression in human breast cancer cells: insights into gene networks and pathways underlying estrogenic control of proliferation and cell phenotype. Endocrinology (2003) 144(10):4562–74. doi: 10.1210/en.2003-0567

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Rothenberger NJ, Somasundaram A, Stabile LP. The role of the estrogen pathway in the tumor microenvironment. Int J Mol Sci (2018) 19(2). doi: 10.3390/ijms19020611

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Weir HK, Marrett LD, Cokkinides V, Barnholtz-Sloan J, Patel P, Tai E, et al. Melanoma in adolescents and young adults (ages 15-39 years): united states, 1999-2006. J Am Acad Dermatol (2011) 65(5 Suppl 1):S38–49. doi: 10.1016/j.jaad.2011.04.038

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Enninga EAL, Moser JC, Weaver AL, Markovic SN, Brewer JD, Leontovich AA, et al. Survival of cutaneous melanoma based on sex, age, and stage in the united states, 1992-2011. Cancer Med (2017) 6(10):2203–12. doi: 10.1002/cam4.1152

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Donley GM, Liu WT, Pfeiffer RM, McDonald EC, Peters KO, Tucker MA, et al. Reproductive factors, exogenous hormone use and incidence of melanoma among women in the united states. Br J Cancer. (2019) 120(7):754–60. doi: 10.1038/s41416-019-0411-z

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Svoronos N, Perales-Puchalt A, Allegrezza MJ, Rutkowski MR, Payne KK, Tesone AJ, et al. Tumor cell-independent estrogen signaling drives disease progression through mobilization of myeloid-derived suppressor cells. Cancer Discovery (2017) 7(1):72–85. doi: 10.1158/2159-8290.CD-16-0502

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Marzagalli M, Montagnani Marelli M, Casati L, Fontana F, Moretti RM, Limonta P. Estrogen receptor β in melanoma: from molecular insights to potential clinical utility. Front Endocrinol (Lausanne). (2016) 7:140. doi: 10.3389/fendo.2016.00140

PubMed Abstract | CrossRef Full Text | Google Scholar

35. de Giorgi V, Mavilia C, Massi D, Gozzini A, Aragona P, Tanini A, et al. Estrogen receptor expression in cutaneous melanoma: a real-time reverse transcriptase-polymerase chain reaction and immunohistochemical study. Arch Dermatol (2009) 145(1):30–6. doi: 10.1001/archdermatol.2008.537

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ning WR, Jiang D, Liu XC, Huang YF, Peng ZP, Jiang ZZ, et al. Carbonic anhydrase XII mediates the survival and prometastatic functions of macrophages in human hepatocellular carcinoma. J Clin Invest (2022) 132(7). doi: 10.1172/JCI153110

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Guo X, Yu X, Li F, Xia Q, Ren H, Chen Z, et al. Identification of survival-related metabolic genes and a novel gene signature predicting the overall survival for patients with uveal melanoma. Ophthalmic Res (2022) 65(5):516–28. doi: 10.1159/000524505

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Tworkoski K, Singhal G, Szpakowski S, Zito CI, Bacchiocchi A, Muthusamy V, et al. Phosphoproteomic screen identifies potential therapeutic targets in melanoma. Mol Cancer Res (2011) 9(6):801–12. doi: 10.1158/1541-7786.MCR-10-0512

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Xu T, Xu W, Zheng Y, Li X, Cai H, Xu Z, et al. Comprehensive FGFR3 alteration-related transcriptomic characterization is involved in immune infiltration and correlated with prognosis and immunotherapy response of bladder cancer. Front Immunol (2022) 13:931906. doi: 10.3389/fimmu.2022.931906

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Filippou PS, Ren AH, Soosaipillai A, Safar R, Prassas I, Diamandis EP, et al. Kallikrein-related peptidases protein expression in lymphoid tissues suggests potential implications in immune response. Clin Biochem (2020) 77:41–7. doi: 10.1016/j.clinbiochem.2019.12.015

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Panos M, Christophi GP, Rodriguez M, Scarisbrick IA. Differential expression of multiple kallikreins in a viral model of multiple sclerosis points to unique roles in the innate and adaptive immune response. Biol Chem (2014) 395(9):1063–73. doi: 10.1515/hsz-2014-0141

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Salerno EP, Bedognetti D, Mauldin IS, Deacon DH, Shea SM, Pinczewski J, et al. Human melanomas and ovarian cancers overexpressing mechanical barrier molecule genes lack immune signatures and have increased patient mortality risk. Oncoimmunology (2016) 5(12):e1240857. doi: 10.1080/2162402X.2016.1240857

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ware AW, Harris JJ, Slatter TL, Cunliffe HE, McDonald FJ. The epithelial sodium channel has a role in breast cancer cell proliferation. Breast Cancer Res Tr. (2021) 187(1):31–43. doi: 10.1007/s10549-021-06133-7

CrossRef Full Text | Google Scholar

44. Lou JY, Wei LJ, Wang H. SCNN1A overexpression correlates with poor prognosis and immune infiltrates in ovarian cancer. Int J Gen Med (2022) 15:1743–63. doi: 10.2147/IJGM.S351976

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Dou YL, Xing JJ, Kong GC, Wang GC, Lou XH, Xiao X, et al. Identification of the E3 ligase TRIM29 as a critical checkpoint regulator of NK cell functions. J Immunol (2019) 203(4):873–80. doi: 10.4049/jimmunol.1900171

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Han J, Zuo J, Zhang X, Wang L, Li D, Wang YD, et al. TRIM29 is differentially expressed in colorectal cancers of different primary locations and affects survival by regulating tumor immunity based on retrospective study and bioinformatics analysis. J Gastrointest Oncol (2022) 13(3):1132–51. doi: 10.21037/jgo-22-365

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Nalbandian G, Paharkova-Vatchkova V, Mao A, Nale S, Kovats S. The selective estrogen receptor modulators, tamoxifen and raloxifene, impair dendritic cell differentiation and activation. J Immunol (2005) 175(4):2666–75. doi: 10.4049/jimmunol.175.4.2666

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Liang X, Briaux A, Becette V, Benoist C, Boulai A, Chemlali W, et al. Molecular profiling of hormone receptor-positive, HER2-negative breast cancers from patients treated with neoadjuvant endocrine therapy in the CARMINA 02 trial (UCBG-0609). J Hematol Oncol (2018) 11(1):124. doi: 10.1186/s13045-018-0670-9

PubMed Abstract | CrossRef Full Text | Google Scholar