Abstract
Melanoma is one of the major cancer types for which new immune-based cancer treatments have achieved promising results. However, anti–PD-1 and anti–CTLA-4 therapies are effective only in some patients. Hence, predictive molecular markers for the development of clinical strategies targeting immune checkpoints are needed. Using The Cancer Genome Atlas (TCGA) RNAseq data, we found that expression of ESRP1, encoding a master splicing regulator in the epithelial–mesenchymal transition (EMT), was inversely correlated with tumor-associated immune cytolytic activity. That association holds up across multiple TCGA tumor types, suggesting a link between tumor EMT status and infiltrating lymphocyte activity. In melanoma, ESRP1 mainly exists in a melanocyte-specific truncated form transcribed from exon 13. This was validated by analyzing CCLE cell line data, public CAGE data, and RT-PCR in primary cultured melanoma cell lines. Based on ESRP1 expression, we divided TCGA melanoma cases into ESRP1-low, -truncated, and –full-length groups. ESRP1-truncated tumors comprise approximately two thirds of melanoma samples and reside in an apparent transitional state between epithelial and mesenchymal phenotypes. ESRP1 full-length tumors express epithelial markers and constitute about 5% of melanoma samples. In contrast, ESRP1-low tumors express mesenchymal markers and are high in immune cytolytic activity as well as PD-L2 and CTLA-4 expression. Those tumors are associated with better patient survival. Results from our study suggest a path toward the use of ESRP1 and other EMT markers as informative biomarkers for immunotherapy. Cancer Immunol Res; 4(6); 552–61. ©2016 AACR.
Introduction
Melanoma is one of the most invasive malignancies, and patients with melanoma mainly die from tumor metastasis. The reversible phenotype switch between proliferative and invasive cells that drives the metastasis of tumor cells has been proposed as a model for melanoma progression (1). Many transcription factors and signaling pathways induced by changes in the microenvironment have been implicated in tumorigenesis and this phenotypic switch (2, 3). MITF, a master transcriptional regulator involved in melanocyte differentiation and pigment production, has been shown to play a critical role in melanoma development (4, 5). High expression of MITF is associated with transformation of melanocytes and hyperproliferation of the transformed cells. Downregulation of MITF expression leads to greater invasive potential of melanoma cells (3, 6, 7). Changes in the Wnt signaling pathway have also been implicated in the phenotypic switch (8–11). The inversely expressed genes ROR1 and ROR2 are reported to be associated with proliferative and invasive phenotypes in melanoma, respectively (12). However, many of these observations are based on in vitro gain-of-function or loss-of-function experiments and might not reflect the complexity of microenvironment changes and cohesive signaling networks in vivo.
ESRP1 encodes a cell type–specific splicing regulator protein exclusively expressed in epithelial cells such as those comprised of ectoderm-derived tissues and endothelial organ lining. Expression of ESRP1 controls cell type–specific alternative splicing of many genes, such as FGFR2, CD44, and CTNND1, that are important signaling molecules mediating cell division, growth, and differentiation under specific microenvironment settings (13). For example, two alternative forms of FGFR2, IIIb and IIIc, which exhibit different affinity with their FGF ligands, are predominantly expressed in epithelial and mesenchymal cells, respectively. The selective expression of the FGFR2 IIIb or IIIc isoforms is under the control of ESRP1 (14). We previously reported that dominance of epithelial or mesenchymal cell types exists in renal cell carcinoma (RCC; ref. 15). Our results support the notion that clear cell renal cell carcinoma (ccRCC) originates from mesenchymal-like stem cells embedded in adult kidney, unlike other subtypes of RCC, such as chromophobe, that likely have an epithelial-like stem cell origin. Thus, ESRP1 and status of its target genes might serve as molecular markers of tumor cells with epithelial or mesenchymal properties.
Cutaneous melanoma originates from melanocytes derived from multipotent neural crest cells (NCC). NCCs undergo epithelial–mesenchymal transition (EMT) before delaminating from the dorsal neural tube and can differentiate into many cell types, such as neurons of the peripheral nervous system, glial and Schwann cells, melanocytes, and craniofacial cartilage. Although the EMT process enables NCC migration, formation of NCCs' derivative structures requires re-aggregation of cells, often coupled with the reverse mesenchymal–epithelial transition (MET) following migration (16, 17). Behavior of neural crest stem cells during embryogenesis and their existence in adult NCC-derived tissues raise questions about the origin and ultimate fate of melanoma cells (18).
Transformation of melanocytes is driven by somatic changes in cancer genomes. However, these somatic changes are not sufficient to fully explain the highly heterogeneous and metastatic nature of the disease. A comprehensive genomic landscape of alterations in cutaneous melanoma was developed based on multidimensional high-throughput data (19). The substantial heterogeneity of tumor cells and their microenvironments is reflected in the complex patterns of genomic, epigenomic, and transcriptomic alternations. With rapid advances of immunotherapy in melanoma and other types of cancers, we urgently need to discover diagnostic and prognostic markers that could be applied in clinical settings. In this study, we further utilized datasets generated for melanoma tumors by The Cancer Genome Atlas (TCGA) to investigate ESRP1-based cell-type plasticity in melanoma and its implication in melanoma biology and treatment.
Materials and Methods
Samples and datasets
TCGA RNA-Seq expression and methylation (HumanMethylation450 probes) data and corresponding clinical metadata were accessed from TCGA public access data portal released before January 27, 2015 (https://tcga-data.nci.nih.gov/tcga/dataAccessMatrix.htm). Exon level mapping files generated for Cancer Cell Line Encyclopedia (CCLE) RNA-Seq gene expression were accessed through the ArrayStudio OncoLand portal (http://www.omicsoft.com). Cap analysis gene expression (CAGE) data were queried from FANTOM5 (http://fantom.gsc.riken.jp). Other RNA-Seq data including melanocyte cell lines (GSM1138580 and GSM819489), shMITF knockdown experiments (SRP048577), mouse melanocyte and melanoma samples (SRX478034, SRX478033, ERX386690, and SRX475297) were downloaded from NCBI's Short Reads Archive. Somatic point mutation data were accessed from GDAC Firehose (http://gdac.broadinstitute.org/).
The melanocyte cell line used in real-time PCR was purchased from SicenCell Research Laboratory, Inc., in 2004. Other melanoma cell lines, either commercially available or established from low-passage culture of cutaneous melanoma lesions and authenticated by branches of Ludwig Cancer Research, were acquired between 1998 and 2013. Specifically, Hs294T was acquired from ATCC. T618A, Me246, Me260, and Me275 were established and authenticated at the Ludwig Institute for Cancer Research, Lausanne Branch. The SK-MEL series of melanoma cell lines were established and authenticated at the Ludwig Institute for Cancer Research, New York branch, at the Memorial Sloan Kettering Cancer Center. LM-MEL-62 was established and authenticated at the Ludwig Institute for Cancer Research Melbourne, Austin Branch, at Melbourne, Australia. Cell lines were usually cultured for a maximum of three passages.
Defining ESRP1 groups based on exon expression
Calculation of exon expression is based on exon quantification file; gene expression is based on normalized gene expression values in RSEM (RNA-Seq by Expectation-Maximization). The ESRP1-low cutoff is defined by ESRP1 gene expression less than 125 RSEM after examination of the ESRP1 expression distribution in all samples. Among the rest that have ESRP1 gene expression >125 RSEM, samples with an exon truncation ratio ≥0.65 are defined as ESRP1 truncated; samples with a truncation ratio smaller than 0.65 are grouped as ESRP1 full length. Considering that exons 14 and 15 are alternatively expressed exons, the exon ratio used in assigning ESRP1-based subgrouping is defined as follows: sum (exon 13 + exon 16) RSEM/sum (exon (1–12) + exon 13 + exon 16) RSEM. The truncation ratio cutoff of 0.65 was defined after examination of its distribution. The heatmap in Fig. 1 was generated without median removal. Median removal was used in all other expression heatmaps.
ESRP1 expression is inversely correlated with immune cytolytic activity in TCGA carcinomas. Cytolytic activity measurement was undertaken as described in Materials and Methods. x-axis, ESRP1 expression in log2 scale; y-axis, cytolytic activity in log2 scale. Primary and metastatic tumors have been included but samples that do not express ESRP1 (RSEM < 1) have been removed. BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; COAD, colon adenocarcinoma; LUSC, lung squamous carcinoma; PAAD, pancreatic adenocarcinoma; THCA, thyroid carcinoma; TNBC, triple-negative breast cancer.
ESRP1 expression is inversely correlated with immune cytolytic activity in TCGA carcinomas. Cytolytic activity measurement was undertaken as described in Materials and Methods. x-axis, ESRP1 expression in log2 scale; y-axis, cytolytic activity in log2 scale. Primary and metastatic tumors have been included but samples that do not express ESRP1 (RSEM < 1) have been removed. BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; COAD, colon adenocarcinoma; LUSC, lung squamous carcinoma; PAAD, pancreatic adenocarcinoma; THCA, thyroid carcinoma; TNBC, triple-negative breast cancer.
Gene expression signatures
Cytolytic activity is calculated as the geometric mean of GZMA and PRF1 using log2 expression values. We used normal tissue gene expression data obtained from the GTEx portal (http://www.gtexportal.org/home/) to identify two skin-specific genes (KRT2 and LOR). With the highest KRT2/LOR expressing non-skin tissue, its median expression of KRT2/LOR plus 3× standard deviation is lower than RPKM (reads per kilobase of transcript per million mapped reads) 16. We therefore call melanoma samples with KRT2 or LOR expression above RPKM 16 as keratinocyte contaminated.
Quantitative real-time reverse transcription PCR
Analysis of cDNA samples was performed in duplicate for ESRP1 and for the reference gene within each experiment using the CFX96 Real-Time PCR system (Biorad) and Taqman platform (Applied Biosystems). TFRC was amplified as an internal reference gene. The PCR primers and probes for ESRP1 (Hs00214472_m1, flanking exons 11–2; Hs00936411_m1, flanking exons 13–14) and TFRC internal control gene (4326323E) were purchased from Applied Biosystems. Primers used for PCR amplification were chosen to encompass intron between exon sequences to avoid amplification of genomic DNA.
Cloning of full-length ESRP1 and transfection into SKMEL28
PCR was undertaken with High Fidelity Taq polymerase (Invitrogen) using cDNA from the BT20 cell line: Forward 5′-ATGACGGCCTCTCCGGATTACTTGGTG-3′ and reverse 5′-AATACAAACCCATTCTTTGGGT-3′. The PCR product was recovered and cloned in pCDNA pcDNA3.1/V5-His using the pcDNA3.1/V5-His TOPO TA Expression Kit (Invitrogen). The correct full-length clone was determined by Sanger sequencing using the T7 promoter and V5 reverse primer.
Transfection was performed using Effectene (Qiagen) according to the manufacturer's instructions. Western blotting was performed as previously described (15). The antibodies used included: mouse monoclonal anti-V5 (clone 2F11F7; Life Technologies) and anti–ESRP-1/2 (clone 23A7.C9; Rockland Immunochemicals) and rabbit polyclonal anti-actin (20-33; Sigma-Aldrich).
Statistical tests
The Krustal–Wallis test (nonparametric ANOVA) was used for multigroup significance tests. Breslow thickness was categorized in four levels as I, < 1 mm; II, < 2 mm; III, < 4 mm; and IV, 4+ mm. The Fisher exact test is used for testing sample distribution. The Kaplan–Meier log-ranked test was performed for survival benefits. All tests were done in R, Graphpad, SSPSv12 or ArrayStudio.
Results
Inverse correlation of ESRP1 expression to tumor-associated immune cytolytic activity
The use of antibodies to CTLA-4 or PD-1 has induced remarkable and durable antitumor responses in patients with melanoma and other types of cancers (20). However, only a minority of patients treated achieve complete tumor regression (21, 22). Therefore, identifying predictive immunotherapy biomarkers is critical for both selecting appropriate patients for immunotherapy and for our complete understanding of the mechanism of action of these agents. It has been reported that PD-L1 expression is regulated by microRNA-200/ZEB1 in lung cancer (23), suggesting that tumor EMT status may be associated with immunotherapy response. To examine the relationship between EMT and tumor immune cytolytic activity, we interrogated multiple cancer types in the TCGA database.
We used ESRP1 as a biomarker of EMT status because its expression is highly restricted to epithelial cells and its level is tightly linked to the selective expression of epithelial or mesenchymal-specific alternative splice forms of multiple genes such as FGFR2 (15, 24). We also applied a two-gene expression signature (PRF1 and GZMA) of immune cytolytic activity (25). These genes encode perforin and granzyme A proteins, which are specifically expressed by CD8+ cytotoxic T cells upon activation. The two-gene signature significantly correlates with T-cell markers, including cytotoxic T lymphocyte (CTL) markers, as well as with pan-cancer survival benefits. As shown in Fig. 1, ESRP1 expression was inversely correlated to this signature of cytolytic activity in many cancer types, including lung, colon, prostate, breast, bladder, and thyroid, with significant P values (< 1 × 10−10; Supplementary Table S1). Although a trend of inverse correlation exists in other cancer types, the Pearson correlation coefficient has shown less statistical significance (Supplementary Table S1).
Tissue-specific expression of truncated ESRP1 transcript in melanoma and normal melanocyte
Intriguingly, melanoma, a tumor of neural crest origin (26), had high ESRP1 expression compared with other neural crest–derived tumors, such as glioblastoma and low-grade glioma (Supplementary Fig. S1A). This was also observed in cancer cell lines [Cancer Cell Line Encyclopedia (CCLE) dataset], where ESRP1 was highly expressed in cells of epithelial origin, as is the case in most carcinoma cell lines, but not expressed in cells of mesenchymal origin, such as sarcomas and lymphomas, nor in other neuronal tumor cell lines (Supplementary Fig. S1B). ESRP1 expression was high in cell lines of “skin” melanoma origin (61/62 CCLE skin cell lines are from melanoma).
A close examination of ESRP1 expression at the exon level in melanoma revealed that a majority of tumors do not express full-length ESRP1 transcript. Instead, only the last four exons (exons 13–16) were highly expressed and presumably would not generate ESRP1 protein that carries normal function. Based on the absolute ESRP1 expression level and the relative exon 13–16 expression ratio (see Materials and Methods), we classified melanomas in three subgroups (Fig. 2A). In approximately 20% of melanomas, designated ESRP1-low, ESRP1 was expressed at low levels (RSEM value of less than 125). In approximately two thirds of melanomas, the ESRP1 truncated form was prominent (exon ratio above 0.65), whereas about 10% of all tumors express full-length ESRP1 transcripts in substantial amounts (Table 1). Notably, the one normal control sample in the TCGA melanoma dataset fell into the ESRP1-truncated group, suggesting that the truncated transcript form was not tumor specific, but rather produced during normal melanocyte differentiation. We did not observe any aberrant change in the ESRP1 homolog, ESRP2, which also has epithelial-specific expression, although higher ESRP2 expression was associated with the ESRP1 full-length group, as expected (Fig. 2A).
Detection of truncated ESRP1 transcripts in melanoma. A, ESRP1 expression at gene and exon levels in melanoma (tumor samples, n = 471; normal melanocyte, n = 1). Melanoma samples are organized into three groups (see Materials and Methods) and ordered by ESRP1 expression. ESRP1 exon expression heatmap is shown at top, followed by ESRP1 expression values in a barplot. A truncation ratio of 0.65 was used to call full-length ESRP1 tumors (see Materials and Methods). Position of the normal sample is marked by a vertical yellow line. Tumor sites for primary, regional subcutaneous, lymph node, and distant metastasis tumors are marked by dark-blue lines. The “keratin” lane marks potential contamination of keratinocytes (see Materials and Methods). B, expression of ESRP1 at exon level in CCLE cancer cell lines. C, detection of low, truncated, and full-length ESRP1 expression in melanoma cell lines by real-time PCR. Two probe sets amplify regions flanking exon 11–12 and exon 13–14, respectively. Top, relative signal intensity of the two probe sets normalized by normal kidney. Bottom, percentage of exon 13–14 amplicon out of the total signal. HEM, Hermes cell line.
Detection of truncated ESRP1 transcripts in melanoma. A, ESRP1 expression at gene and exon levels in melanoma (tumor samples, n = 471; normal melanocyte, n = 1). Melanoma samples are organized into three groups (see Materials and Methods) and ordered by ESRP1 expression. ESRP1 exon expression heatmap is shown at top, followed by ESRP1 expression values in a barplot. A truncation ratio of 0.65 was used to call full-length ESRP1 tumors (see Materials and Methods). Position of the normal sample is marked by a vertical yellow line. Tumor sites for primary, regional subcutaneous, lymph node, and distant metastasis tumors are marked by dark-blue lines. The “keratin” lane marks potential contamination of keratinocytes (see Materials and Methods). B, expression of ESRP1 at exon level in CCLE cancer cell lines. C, detection of low, truncated, and full-length ESRP1 expression in melanoma cell lines by real-time PCR. Two probe sets amplify regions flanking exon 11–12 and exon 13–14, respectively. Top, relative signal intensity of the two probe sets normalized by normal kidney. Bottom, percentage of exon 13–14 amplicon out of the total signal. HEM, Hermes cell line.
Distribution of melanoma tumors on ESRP1 expression status and tumor site before and after (before/after) removing potentially keratinocyte contaminated samples
. | ESRP1-low . | ESRP1-trun . | ESRP1 full-length . | Total . |
---|---|---|---|---|
Primary tumor | 4a/3a | 69/56 | 32a/5 | 105/64 |
Regional lymph node | 56/52 | 154/153 | 12/12 | 222/217 |
Cutaneous and subcutaneous | 19/19 | 49/47 | 6/3 | 74/69 |
Distant metastasis | 13/13 | 49/47 | 6/5 | 68/65 |
Total | 92 | 321/303 | 56/25 | 469/415 |
. | ESRP1-low . | ESRP1-trun . | ESRP1 full-length . | Total . |
---|---|---|---|---|
Primary tumor | 4a/3a | 69/56 | 32a/5 | 105/64 |
Regional lymph node | 56/52 | 154/153 | 12/12 | 222/217 |
Cutaneous and subcutaneous | 19/19 | 49/47 | 6/3 | 74/69 |
Distant metastasis | 13/13 | 49/47 | 6/5 | 68/65 |
Total | 92 | 321/303 | 56/25 | 469/415 |
a, P value < 0.0001 by the Fisher exact test.
When tumor sites were aligned to ESRP1 groups, we observed enrichment of primary tumors inside the ESRP1 full-length group (Fig. 2A and Table 1). This raised a concern about whether full-length ESRP1 expression was from contaminating keratinocytes, which are often found in primary melanomas. Using gene expression data from normal tissue, obtained from the GTEx portal (http://www.gtexportal.org/home/), we identified two skin-specific genes (KRT2 and LOR) with strong expression (Supplementary Fig. S2A), and used these two genes to estimate keratinocyte contamination in TCGA samples. Indeed, about half of ESRP1 full-length tumors were KRT2 and/or LOR positive. After removing these contaminated samples, a portion of melanomas still expressed full-length ESRP1. These melanomas had been collected from different tumor sites, including primary, regional subcutaneous, lymph node, and distant metastasis sites, with no enrichment of any particular tumor site (Supplementary Fig. S2B). Similar findings were observed for ESRP1 expression in the CCLE cell line database, in which 4 out of 53 melanoma cell lines had full-length ESRP1 expression and 10 cell lines had low ESRP1 expression (Fig. 2B). Therefore, although ESRP1 full-length expression in some samples may be introduced by keratinocyte contamination, some melanoma tumors did express full-length ESRP1. To avoid signals from contaminating keratinocytes, we discarded KRT2/LOR-positive samples from further analysis (Table 1).
We further validated our findings by real-time PCR with a separate collection of 42 melanoma cell lines derived from cultures of primary tumors and a melanocyte cell line (see Materials and Methods). To detect the truncated transcripts, we used two sets of PCR primers to amplify cDNA flanking ESRP1 exons 11 and 12 or exons 13 and 14, respectively, and determine the relative abundance of exon 13–14 versus exon 11–12 products. In comparison with normal kidney and colon samples, all melanoma cell lines except one (SKMEL128), as well as the normal melanocyte cell line (HEM), exhibited an increased usage of exons 13–14, confirming that the truncated ESRP1 transcript predominates in the majority of melanomas (Fig. 2C). Because in the TCGA dataset we only detected expression of this truncated form in melanoma, and the presence of this particular transcript in other types of cancers is rare (Supplementary Fig. S1B), expression of the truncated ESRP1 was apparently melanocyte specific.
Mechanism of truncated ESRP1 expression in melanoma
To further investigate the truncated form of ESRP1, we performed additional analyses. First, we tested its origin by cloning full-length cDNA of ESRP1 into an expression vector and transiently transfecting into SKMEL28, a melanoma cell line that predominantly expresses the truncated ESRP1 transcripts (Fig. 2C). Only full-length ESRP1 protein was detected in the transfected cell lysate by Western blotting after 48 and 72 hours (data not shown). This suggests that SKMEL28 cells could express and sustain expression of full-length ESRP1 and the truncated transcript was not the post-transcriptional processed product from the full-length ESRP1. Because the region of ESRP1 recognized by the monoclonal antibody raised against the full-length mouse Esrp1 is unknown, we could not exclude the possibility that the truncated ESRP1 protein was also produced in these cells.
We further hypothesized that an alternative transcriptional start site (TSS) could direct expression of the truncated ESRP1 transcript. To test this assumption, we queried high-throughput sequencing mapping outputs of CAGE libraries, which capture the 5′ end sequence of mRNAs in the FANTOM project (http://fantom.gsc.riken.jp). Indeed, in two skin melanoma cell lines (COLO-679 and G-361), most CAGE reads for ESRP1 were mapped to a region inside exon 13 (Fig. 3A). As a control, CAGE reads were seen to be aligned at the 5′ end of MLANA (the melanoma-specific melan-A gene) promoter region, as expected. In other epithelial cell lines, CAGE reads were almost exclusively mapped to ESRP1 5′ end promoter region. Thus, the truncated ESRP1 resulted from alternative TSS usage at an approximate location of chr8:95690440 (hg19) at the start of exon 13 in melanoma. The two known functional domains (RNaseH and RRM domains) are missing from the truncated ESRP1, so the product from this new isoform is unlikely to retain its normal function (Fig. 3B).
Melanomas express a novel transcriptional isoform of ESRP1. A, alternative TSS of ESRP1. Top, CAGE libraries capture reads mapped to TSS sites at MLANA (top) and ESRP1 (bottom) gene loci. Depth of reads coverage is labeled on the right. Exons in alternative transcripts are drawn in filled boxes. B, ESRP1 coding structure. The E-box location is marked out by a red box. Exon 13 is labeled to show the relative exon positions. Alternative TSS for the truncated ESRP1 within exon 13 is marked by a red arrow. Functional domains encoded by ESRP1, including RNaseH-like domain and RRM (RNA recognition motif) are shown in colored boxes in the gene structure. Sequence surrounding the junction of intron 12 and exon 13 is shown below with comparison among other species. C, correlation between MITF and ESRP1 gene expression. Samples are categorized into ESRP1-low and ESRP1-trun groups. D, shMITF downregulates expression of the truncated ESRP1. CDK2, a known MITF target, and ACTB (Actin beta) are used as positive and negative controls of the experiment, respectively.
Melanomas express a novel transcriptional isoform of ESRP1. A, alternative TSS of ESRP1. Top, CAGE libraries capture reads mapped to TSS sites at MLANA (top) and ESRP1 (bottom) gene loci. Depth of reads coverage is labeled on the right. Exons in alternative transcripts are drawn in filled boxes. B, ESRP1 coding structure. The E-box location is marked out by a red box. Exon 13 is labeled to show the relative exon positions. Alternative TSS for the truncated ESRP1 within exon 13 is marked by a red arrow. Functional domains encoded by ESRP1, including RNaseH-like domain and RRM (RNA recognition motif) are shown in colored boxes in the gene structure. Sequence surrounding the junction of intron 12 and exon 13 is shown below with comparison among other species. C, correlation between MITF and ESRP1 gene expression. Samples are categorized into ESRP1-low and ESRP1-trun groups. D, shMITF downregulates expression of the truncated ESRP1. CDK2, a known MITF target, and ACTB (Actin beta) are used as positive and negative controls of the experiment, respectively.
Upstream sequence analysis identified a consensus E-box motif (CACGTG) located −57-bp upstream of ESRP1 exon 13 (Fig. 3B), which is a known binding site for MITF, a master regulator of melanocyte development and a known melanoma oncogene (4, 5). Several lines of evidence support the idea that MITF regulated the expression of truncated ESRP1 through binding to the E-box motifs inside intron 12 and initiated new transcription from within exon 13. First, correlation of global gene expression to truncated ESRP1 expression resulted in MITF as the most correlated gene, with a Pearson coefficient of 0.77 (Fig. 3C), and with many of the remaining top correlated genes being MITF target genes. Second, the consensus E-box found in humans is conserved in monkeys and apes but not in mice and rats (Fig. 3B and Supplementary Fig. S3). Consistently, we detected little mouse Esrp1 expression using mouse nevi and melanoma RNA-Seq samples (see Materials and Methods). Third, knockdown of MITF by shRNA in 501MEL and Hermes cell lines leads to reduced expression of MITF target genes (e.g. CDK2) and of truncated ESRP1 (Fig. 3D). Together, these results suggest that MITF regulates truncated ESRP1 expression in melanomas.
Epithelial and mesenchymal cell phenotypes are associated with ESRP1 status in melanoma
We examined the expression of a set of epithelial and mesenchymal markers in the three ESRP1 subgroups after removing keratinocyte-contaminated samples (Fig. 4A). ESRP1-low tumors exhibited a strong expression signature of mesenchymal markers like N-cadherin (CDH2) and TWIST1. About half of melanomas with full-length ESRP1 expression displayed higher expression of an epithelial marker panel including genes such as for keratin. This was not restricted to primary tumors, but was also found in tumors from regional cutaneous/subcutaneous locations, lymph nodes, and distant metastatic sites. Approximately two thirds of tumors expressing truncated ESRP1 seemed to reside in a state between epithelial and mesenchymal, even though ESRP1 was functionally compromised in these tumors. This suggested that loss of ESRP1 function per se is not sufficient to render a mesenchymal phenotype, at least in melanoma. Surprisingly, E-cadherin (CDH1) had higher expression in both ESRP1 full-length– and ESRP1-truncated groups. We also examined the isoform switch in the FGFR2 gene, which showed a clear transition from the epithelial isoform (FGFR2-b) in ESRP1 full-length tumors to the mesenchymal isoform (FGFR2-c) in tumors with truncated and low ESRP1 (Fig. 4B). This finding not only validated our computational prediction of ESRP1 full-length expression, but also confirmed that truncated ESRP1 was unable to exert its normal function to induce isoform switch in the FGFR2 gene.
ESRP1 subgroups are associated with epithelial–mesenchymal cellular properties. A, expression heatmap of molecular markers for epithelial and mesenchymal cells. B, FGFR2 expression at exon level shows that the mesenchymal and epithelial forms are dominantly expressed in non–ESRP1 full-length and ESRP1 full-length subgroups, respectively. C, heatmap of the 105-gene signature and MITF expression levels (RSEM) are aligned to ESRP1 subgroups. Samples in all heatmaps are in the same order as in Fig. 2.
ESRP1 subgroups are associated with epithelial–mesenchymal cellular properties. A, expression heatmap of molecular markers for epithelial and mesenchymal cells. B, FGFR2 expression at exon level shows that the mesenchymal and epithelial forms are dominantly expressed in non–ESRP1 full-length and ESRP1 full-length subgroups, respectively. C, heatmap of the 105-gene signature and MITF expression levels (RSEM) are aligned to ESRP1 subgroups. Samples in all heatmaps are in the same order as in Fig. 2.
Previously, several studies defined molecular subtypes of melanoma using gene expression profiling (3, 27–29). Among those, Hoek and colleagues classified melanoma into “proliferative” and “invasive” subtypes based on the microarray profiling of 86 cultured melanomas (3). When we compared the 105-gene signature from that study with our ESRP1 subgroups, we observed a good correlation of the ESRP1-low group to the invasive subtype, as well as the truncated-ESRP1 group to the proliferative subtype (Fig. 4C). The ESRP1 full-length group, however, comprised samples with both proliferative and invasive signatures. MITF is apparently a dominant factor in the proliferative subtype because many MITF target genes were in the 105-gene list and highly expressed in the proliferative subtype, whereas in the ESRP1-low/invasive group, MITF expression was low (Fig. 4C). This implies that increased MITF level may actively block a mesenchymal phenotype, in agreement with previous findings (3).
ESRP1-low melanomas are associated with increased immune cytotoxicity and favorable patient survival
We examined immune cytolytic activity using cytotoxic T-cell marker genes and the two-gene signature in three ESRP1-based melanoma subgroups after removing all primary melanoma samples, because these are not particularly relevant to immune checkpoint therapy. This showed a statistically significant increase of cytolytic activity in the ESRP1-low melanomas (Fig. 5A), which was consistent with our findings in other cancer types (Fig. 1). Predicted cytolytic activity between the ESRP1 full-length and the ESRP1-truncated groups was not significantly different, suggesting the correlation of cytolytic activity was more related to the mesenchymal phenotype rather than to the epithelial phenotype. We next examined the expression of the PD-1, PD-L1, PD-L2, and CTLA-4 genes in ESRP1 subgroups (Fig. 5A and B). All these molecules, which are targeted by immune checkpoint blockage antibodies, had increased expression in the ESRP1-low tumors. We also examined the possible relationship of ESRP1 grouping to mutational status and clinical parameters. No particular association was found between BRAF and NRAS mutation and ESRP1 groups, although KIT mutations seem to be absent in ESRP1-low tumors (Fig. 5A). The ESRP1 groups were also not associated with the total number of somatic mutations. In the TCGA melanoma dataset, the ESRP1 grouping had a moderate association to Breslow depth, with ESRP1-low tumors having the thinnest Breslow depth (Fig. 5B). Kaplan–Meier survival analysis on TCGA melanoma patients stratified by ESRP1 status showed a trend of more favorable survival of ESRP1-low tumors (P = 0.057, Fig. 5C), which is consistent with the lower Breslow depth values and higher active immune signatures found within this group.
Association of ESRP1 expression status with tumor-associated cytolytic activity and clinical parameters. A, comparison of immune cytolytic activity and related genes in ESRP1 subgroups. Top, heatmap of expression of immune checkpoint target genes and active T-cell markers. Bottom, alignment of mutations and clinical parameters to ESRP1 subgroups. B, boxplot presentation of expression levels of PD-L2, CTLA-4, cytolytic activity, Breslow depth, and tumor T stages in ESRP1 subgroups. C, Kaplan–Meier survival analysis of TCGA melanoma cases stratified by ESRP1 status.
Association of ESRP1 expression status with tumor-associated cytolytic activity and clinical parameters. A, comparison of immune cytolytic activity and related genes in ESRP1 subgroups. Top, heatmap of expression of immune checkpoint target genes and active T-cell markers. Bottom, alignment of mutations and clinical parameters to ESRP1 subgroups. B, boxplot presentation of expression levels of PD-L2, CTLA-4, cytolytic activity, Breslow depth, and tumor T stages in ESRP1 subgroups. C, Kaplan–Meier survival analysis of TCGA melanoma cases stratified by ESRP1 status.
In summary, our data show an inverse correlation between ESPR1 expression and tumor-associated immune cytolytic activity in multiple tumor types. Through the study of altered ESRP1 transcription in melanoma, we gained further insights on the complex heterogeneity of melanoma cell properties during tumor progression, and strengthened a link between tumor-intrinsic EMT status and tumor microenvironment in melanoma. Results from our study highlight the potential utility of ESRP1 status in predicting response to checkpoint blockade immunotherapy. The conceptual link between tumor EMT status and activated infiltrating lymphocytes may have further implications in the immunotherapy field and will need further elucidation.
Discussion
Cancer immunotherapy by immune checkpoint inhibitors, adoptive T-cell therapy, and therapeutic vaccines has advanced to the first line of treatments for many metastatic cancers, especially melanoma (20, 30, 31). Contributing to this opportunity is the rich stromal cell–orchestrated tumor microenvironment, particularly those tumors with high numbers of tumor-infiltrating lymphocytes (TIL), as is the case in melanoma. Here, we show that melanoma mesenchymal-like tumors were associated with enhanced immune cytolytic activity. Moreover, this phenomenon was corroborated by a negative correlation between immune cytolytic activity and epithelial phenotypes, and especially with the fact that ESRP1 expression was inversely associated with higher cytolytic activity in many other cancer types (Fig. 1). RCC is another cancer type that has shown relatively high sensitivity to immunotherapy (32). Previously, we have demonstrated that ccRCC has a more mesenchymal-like phenotype (15). Our current results suggest that tumor cells of mesenchymal nature are able to better acquire TILs and are well-suited to immunotherapy. With the success of immune checkpoint blockage therapy (anti–PD1/PDL-1 and/or anti–CTLA-4) in some patients, an important issue is why some patients are either not responsive or later develop resistance. Among the many factors that contribute to tumor immunity, such as HLA abundance, neoantigen load, and intrinsic properties, our findings suggest that the cellular state also has an impact on immune response susceptibility.
In this study, we have described the identification of a melanocyte-specific ESRP1 mRNA transcript. This highly expressed isoform is a result of de novo transcription starting from exon 13, leading to a shorter mRNA (truncated ESRP1) that comprises only the last four exons and loses the region coding for the two functional domains of the protein. Although the TCGA melanoma dataset contains only one normal control sample, which showed truncated ESRP1 expression, we did confirm this finding by RT-PCR in a normal melanocyte cell line (HEM; Fig. 2C). Additionally, publicly accessible RNA-Seq data of independently derived melanocyte specimens (NCBI GEO GSE46805 and GSE33092) all show truncated ESRP1 expression in melanocytes (33, 34).
The TCGA melanoma dataset is a heterogeneous collection of samples from primary tumors, regional, lymph nodes, and distant metastases. Our initial findings that primary tumors were significantly enriched in ESRP1 full-length tumors (Fig. 2A) raised the concern of keratinocyte contamination as the cause of ESRP1 full-length expression. Using a two-gene signature (KRT2/LOR), we were able to predict and remove a majority if not all keratinocyte contaminated samples from our study. KRT2 encodes keratin 2E, an epidermal keratin expressed largely in the upper spinous layer of epidermal keratinocytes. LOR encodes loricrin, which is a major protein component of the cornified cell envelope found in terminally differentiated epidermal cells. We show that expression of KRT2 and LOR is correlated and highly skin specific (Supplementary Fig. S2A), and the function of these genes suggests they are authentic keratinocyte genes. Even after removal of contaminated samples, a small portion of melanomas still expressed full-length ESRP1, which suggests that expression of full-length ESRP1 could be an intrinsic tumor property.
Based on ESRP1 status, we divided melanomas into three subgroups. Approximately two thirds of melanomas expressed the truncated ESRP1 form that was also observed in normal controls, so this likely represents the default ESRP1 status. Based on the expression of EMT markers, the state of these samples is neither distinctly mesenchymal nor epithelial. Compared with loss of expression of other epithelial genes, the E-cadherin (CDH1) gene is still highly expressed in this subgroup. We have observed that some tumors lose ESRP1 expression and are mesenchymal, and we show that this may be linked to the loss of MITF expression. For other tumors, which have epithelial characteristics, full-length ESRP1 is expressed. The underlying basis and mechanisms for these alterations will need further investigation. However, our preliminary results revealed loss of ESRP1 promoter methylation in the ESRP1 full-length group, indicating the involvement of epigenetic regulation. Also, a few samples that expressed full-length ESRP1 did not highly express epithelial markers. A closer examination of these individual samples identified at least one sample that expressed another novel truncated form, starting from exon 3 of ESRP1, and therefore had escaped our initial computational screen. Consequently, other factors may influence the functional status of ESRP1. The ESRP1-low subgroup was associated with better TILs and overall survival compared with the other two subgroups.
The identification of effective markers for diagnosis and prediction of treatment response in cancer immunotherapy is important to help guide effective intervention strategies for individual patients. Our study indicates that ESRP1-low melanomas with greater immune cytolytic activity have a higher probability of responding to immune checkpoint blockade therapy and suggests a path toward the use of ESRP1 and other EMT markers as informative biomarkers for immunotherapy.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: J. Yao, R.L. Strausberg, Q. Zhao
Development of methodology: J. Yao, O.L. Caballero, J.S. Cebon, Q. Zhao
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): O.L. Caballero, D. Rimoldi, A. Behren, J.S. Cebon, Q. Zhao
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Yao, Y. Huang, J.S. Cebon, Q. Zhao
Writing, review, and/or revision of the manuscript: J. Yao, O.L. Caballero, C. Lin, A. Behren, J.S. Cebon, R.L. Strausberg, Q. Zhao
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Behren, M.-C. Hung, J.N. Weinstein
Study supervision: J.S. Cebon, Q. Zhao
Grant Support
This work is funded by Ludwig Cancer Research and Regeneron Pharmaceuticals Inc.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.