Introduction

Endometrial cancer (EMCA) is a clinically heterogeneous disease. Majority of endometrial carcinomas are generally low grade and low stage with favorable prognoses, however, the high-grade EMCA accounts for a disproportionate number of EMCA deaths1,2,3. EMCA has been grouped into 2 types. Type 1 is estrogen potentiated, estrogen receptor (ER) and progesterone receptor (PR) positive, and generally carries a favorable prognosis. Type 2 is ER/PR negative tumors, of non-endometrioid histology (mainly serous and clear cell carcinoma), are seen in post-menopausal women, and are associated with atrophic endometrium, and poor outcomes1. High-grade endometrial carcinoma constitutes a biologically, morphologically, genetically, and clinically heterogeneous group of tumors. Recent developments of large-scale genomic studies reveal that this heterogeneity may be a function of the diversity of various molecular alterations during disease progression. The analyses using The Cancer Genome Atlas (TCGA) data have led to an integrated genomic classification of endometrioid endometrial carcinomas (EECs) and serous endometrial carcinomas (SECs) and the identification of the POLE (ultramutated), microsatellite instability (MSI) (hypermutated), copy-number low (endometrioid) and copy-number high (serous-like) subtypes, with distinct combinations of genomic and epigenetic alterations1. Mounting evidence suggests that some molecular alterations are preferentially found in endometrioid endometrial carcinomas (EECs), including mutations in PTEN and CTNNB1, whereas others such as TP53 mutations are more prevalent in serous endometrial carcinomas (SECs)2,3.

The American Cancer Society estimates that 63,230 new cases of cancer of the body of the uterus (uterine body or corpus) will be diagnosed and about 11,350 women will die from cancers of the uterine body4. Although there are many drugs approved for the treatment of ovarian cancer, there is only one FDA-approved drug (Megestrol Acetate) for EMCA, highlighting the need for new therapies to treat advanced, recurrent and metastatic EMCA5,6. Our laboratory identified nuclear expression of the Yes-associated protein (YAP) as a poor prognostic indicator in the overall survival of patients with EMCA7. YAP, the main downstream target of the Hippo pathway, plays an important role in the balance between cell proliferation and apoptosis8,9,10. Verteporfin (VP)11, an FDA approved drug used in photodynamic therapy (PDT) for adult macular degeneration was recently identified as an inhibitor of YAP and its binding partner TEA Domain Transcription Factor 1 (TEAD) binding12. Since the identification of VP as a YAP/TEAD inhibitor, several in vitro and in vivo studies have revealed the new potential of YAP1 in different cancers, where YAP is overexpressed13,14,15,16. We tested the efficacy of VP treatment in Type 1 EMCA cells (HEC-1-A and HEC-1-B) and observed cytotoxic and anti-proliferative effects17. Based on the molecular heterogeneity observed in EMCA patients and the effects of VP on EMCA cells, we hypothesized that VP might alter the biological processes and pathways associated with progression of EMCA cells. The aims of this study were to study the effects of VP on (1) genes or gene expression modules representative of biological processes known to play a role in EMCA, and (2) to define the association of these genes and/or pathways in the progression of EMCA, using RNA sequencing data. Here, we have used RNA sequencing (RNA-seq) to develop transcriptome data set of control and VP treated EMCA cells. We preferred RNA-seq compared to microarray technologies, as RNA-seq has a better dynamic range in estimates of gene expression and better precision18.

Materials and Methods

EMCA cell lines and culture conditions

We used Type 1 EMCA cells for the RNA-seq analysis portion of this study. HEC-1-A (ATCC, HTB-112) and HEC-1-B (ATCC, HTB-113) were obtained from the American Type Culture Collection (ATCC) (Manassas, VA). Both these cell lines were isolated from a patient with stage IA endometrial cancer. HEC-1-A cells were cultured in McCoy’s 5A medium (ATCC, Manassas, VA) supplemented with 10% (v/v) fetal bovine serum (FBS) (Thermo Fisher Scientific, Waltham, MA), HEC-1-B in Eagle’s minimum essential medium (EMEM) (ATCC, Manassas, VA) supplemented with 10%(v/v) FBS. Antibiotics (10 units/ml of penicillin and 10 mg/ml of streptomycin) were added to all culture media. Both cell lines were incubated at 37 °C in a humidified atmosphere containing 5% carbon dioxide.

Verteporfin (VP) treatment

Verteporfin (Sigma, Cat. No. SML0534) was dissolved in DMSO and added to the medium for a final concentration of 10 nM and the cells were treated for 3 h. Equal concentration of DMSO was added to the control cells.

Sample and library preparation

After VP treatment for 3 h at 10 nM, total RNA was isolated from EMCA cells using RNeasy Plus Mini Kit (Qiagen) with DNase treatment. The RNA concentration and purity were measured using the Nanodrop spectrophotometer (Thermo Fisher Scientific), and RNA integrity was evaluated with the Bioanalyzer RNA 6000 Nano Chip. RNA quality was uniformly excellent and met the following criteria; Nanodrop, 260/280 ratio >1.8; Bioanalyzer, RIN > 6.6. The samples were prepared using Illumina TruSeq Stranded Total RNA Sample Prep Kit with Ribo-Zero Gold and then subjected to 125 cycle paired-end sequencing. The average insert size of libraries constructed with the Illumina TruSeq Stranded Total RNA Sample Prep Kit is approximately 150 bp with inserts ranging from 100 to 400 bp. For the RNA-seq analysis, there were two replicates of the control samples (C) and VP-treated samples (VP) for each cell line (HEC-1-A and HEC-1-B) for a total of 8 samples. Sequencing and analysis were done by High-Throughput Genomics and Bioinformatic Analysis, University of Utah Shared Resources (https://healthcare.utah.edu/huntsmancancerinstitute/research/shared-resources/center-managed/bioinformatics/).

Sequence reads alignment and transcript assembly

After processing on the Illumina HiSeq 2500 instrument, FASTQ files containing the nucleotide scores and quality scores for each position were generated; these reads were aligned to the human reference genome (Ensembl release 75) using the STAR read aligner19, specifying that the sequencing experiment was paired-end. The alignment step produced a SAM file. SAMtools was used to generate BAM files from these for the sequencing runs required to generate count matrices20. After reading in the genomic features model (Ensembl’s GTF file) to count reads/fragments and ensuring none of our sequences were circular, a count of all of the exons grouped by gene was calculated.

Analysis of differentially expressed genes

The Bioconductor RNA-seq workflow was followed to detect differential expression21, using the DESeq 2 and other Bioconductor packages in R22. Because counts for transcripts in RNA-seq data can contain many rows with only zeros, rows with little to no information about the amount of gene expression were removed to reduce the size of the data object and to increase computational speed. When analyzing the general effect of VP treatment in both cell lines, standard RNA-seq workflow was followed in removing 292 features with zero reads and 1087 features with 5 of more reads. For the analysis using only HEC-1-B samples, we removed 656 features with zero counts and 1,776 features with fewer than 5 reads. The correlation plots between the replicates of the cell lines (e.g. 2 replicates of HEC-1-A control) show good correlation of normalized transcript counts (measured in log10(FPKM + 1)). DESeq 2 feature counts and other functions were run on the tables of counts to determine differentially expressed genes before and after VP-treatment. Results were considered statistically significant at an adjusted p < 0.05 (DMSO treated vs VP treated). The table of counts for each condition were used as inputs to the DESeq feature counts function to determine differentially expressed genes before and after VP-treatment. Since the expected variance of RNA-seq counts increases with the mean, DESeq 2’s regularized-logarithm transformation (rlog) of the count data was used to steady the variance across the mean21,22. The regularized log transforms (rlog) were taken of the result to generate a matrix of regularized counts for sample visualizations (Fig. 1). To check gene expression after VP treatment, the top 40 most significant genes were sorted by log fold change (logFC) and clustered in a hierarchical fashion using the rlog differences. This was plotted in volcano plots and heatmaps with accompanying dendrograms using pheatmap and ggplot2 packages in R23. An FDR-adjusted P < 0.05 was considered significant24.

Figure 1
figure 1

Gene expression modulated by VP. (A) Expression heat map of sample-to-sample distances on the matrix of variance-stabilized data for overall gene expression. trt = treatment; C = Control; VP = Verteporfin treated. Clustering differentiates between control and VP treated samples. This heatmap was built using DESeq 2 on normalized gene read counts. All of the rlog values of the dispersion estimates were clustered using the R distance function (dist) to calculate the Euclidean distance between samples. Distance plot for HEC-1A and HEC-1B cell lines showing their control and VP-treated versions. (B) Differential gene expression, with fold difference between log2 normalized expression in control (n = 2) and VP treated (n = 2) plotted versus −log10 adjusted P-value. Each gene is colored based on the log10 base mean expression.

To compare VP-treated and untreated samples, a logFC shrinkage method (apeglm) was used to get true low bias logFC estimates for true large differences25. This shrinkage method uses a Bayesian procedure to moderate fold changes from genes with very low and highly variable counts. This serves to reduce noise in differential transcript counts. The logFC after VP treatment was plotted on the y-axis and the average of the counts normalized by size factor plotted on the x-axis to create an MA-plot with no shrinkage and one with logFC shrinkage (Supplementary Fig. S3).

The top 20 up- and down- regulated genes were also sorted by p-value (adjusted for multiple testing using FDR) and plotted those values in heatmaps as well using the pheatmap package wrapped inside the Huntsman Cancer Institute’s hciR package (Fig. 2). False discovery rates (FDR) were calculated by determining the number of control sample differential expression ratios (i.e. ratios of expression values between replicate controls) that exceeded an VP treated sample/control sample ratio and dividing by the rank order of the VP treated sample/control sample ratio. This is the FDR, i.e. fraction of genes that have higher apparent differential expression ratios by chance. 5% of genes with an FDR ≤ 0.05 are expected to be false discoveries. Only genes for which at least one condition had a differential expression ratio with an FDR ≤ 0.05 were included in further analysis. These were converted to log2 values.

Figure 2
figure 2

Differential RNA-seq feature counts for both EMCA cell lines. (A) Heatmap showing differential RNA-seq feature counts for both EMCA cell lines (HEC-1-A and HEC-1-B) before and after VP treatment sorted by logFC. Heatmap shows genes with RNA expression most altered after VP treatment. (B) Heatmap showing differential RNA-seq feature counts for both HEC-1A and HEC-1B EMCA cell lines combined, before and after VP treatment sorted by adjusted p-value, top 20 up and top 20 down-regulated genes are shown. (C) Differential RNA-seq feature counts for HEC-1-B cells before and after VP treatment sorted by logFC. Heatmap of genes with RNA expression most altered after VP treatment, sorted by log Fold Change (logFC). There were two replicates of the control samples (C) and VP-treated samples (VP) for 4 total HEC-1B samples. (D) Heatmap showing differential RNA-seq feature counts for HEC-1-B cells before and after VP treatment sorted by adjusted p-value, top 20 up and top 20 down-regulated genes.

Protein-protein interaction network

To investigate the interactions between the protein products of top differentially expressed genes in the HEC-1-B cell line following VP treatment, a STRING network was constructed with the nodes consisting of genes and edges derived from experimentally validated protein-protein interactions. From the list of differentially expressed genes, a hypergeometric score was calculated of the likelihood of co-occurrence of genes (Supplementary Table S3). The Pfam protein domains and KEGG pathways most associated with the differentially expressed genes were also obtained through STRING (Supplementary Table S5). Many of the genes altered in HEC-1-B alone seemed to overlap with genes in the cell cycle pathway; a simple overlap test was performed to produce the list of cell cycle pathway genes differentially expressed between VP-treated and Control.

shRNA treatment

The shRNA targeting the YAP gene was used for downregulating YAP. The control siRNA (sc-37007; Santa Cruz) was used as a negative control. shYAP was transfected into EMCA cells (70% confluency) using Lipofectamine3000 (Invitrogen) according to the manufacturer’s instructions. After 72 hours of transfection, the cells were harvested and used for checking the downregulation of YAP on cell cycle proteins. The knockdown of the YAP was verified by western blotting.

Quantitative real time PCR

All primer sequences were determined using established human or mouse GenBank sequences. Primer sequences were designed using PrimerQuest (IDT) software (Supplementary Table S9). For Real-time polymerase chain reaction (RT-PCR) analysis, total RNA was isolated from control and VP (10 nM, 3 h) treated EMCA cells. Total cellular RNA was extracted with RNeasy kit (Qiagen) with on-column DNase treatment. We used RNA whose A260: A280 ratio is ≥2.0. Total RNA was reverse transcribed into first strand cDNA using iScript cDNA Synthesis Kit (Bio-Rad), as per the manufacturer’s instructions. Quantitative analysis of genes was done by SYBR green based real-time PCR using Applied Biosystems Real-Time PCR Detection System. Each sample was measured in triplicate and normalized to the reference GAPDH/β-actin/PGK1/LDHA/PPIH gene expression. A statistical evaluation of RT-PCR results was performed using one-way analysis of variance (ANOVA) to compare test gene expression between control and VP treated cells.

Western blot analysis

Control and VP-treated mice tissues or control and VP-treated cells were lysed in RIPA buffer (Boston Bioproducts, Cat. No. BP-115DG) supplemented with protease and phosphatase inhibitors and subjected to SDS-PAGE. Samples were separated electrophoretically on 10% to 12% gels, electroblotted onto nitrocellulose membrane (Bio-Rad), blots were blocked at room temperature for 1 h in 5% (w/v) milk in phosphate-buffered saline and incubated overnight at 4 °C with primary antibodies. Protein bands were visualized with SuperSignal™ West Pico Chemiluminescent Substrate (Thermo Fisher Scientific) and detected using LAS-3000 (Fujifilm, Tokyo, Japan). All the antibodies used in this study and their details were provided in Supplementary Table S10.

Mice experiments

To evaluate the efficacy of VP for inhibition of tumor growth in a mouse model, we injected HEC-1-B GFP cells into NCr nude mice by subcutaneous (SC) administration. Mice experiments were conducted by Anticancer Inc. (http://www.anticancer.com). All animal studies were conducted with an Anticancer Institutional Animal Care and Use Committee (IACUC)-protocol specifically approved for this study and in accordance with the principles and procedures outlined in the National Institute of Health Guide for the Care and Use of Animals. All animal procedures were carried out under specific-pathogen-free (SPF) conditions. Log-phase HEC-1-B GFP cells (5 × 106) were suspended in 100 µl PBS and injected SC in left flank of mouse. After the tumors reached a size >100 mm3, VP was administered IP at a dose 50 mg/kg bodyweight to VP-treated group mice (n = 10). VP was given 3 times a week for 3 weeks. Control group mice (n = 10) were administered DMSO in a similar manner. Tumor size was measured by caliper and by GFP imaging twice a week. Body weight was measured twice per week and Body condition scores (BCS) were taken daily. Mice were sacrificed after 3 weeks of treatment. Primary tumors were excised and weighed at necropsy. Blood samples were collected, immediately processed to prepare plasma and flash frozen. Based on tumor volume, VP-treated mice were divided into Responders (tumor volume <500 mm3) and Non-responders (tumor volume >500 mm3). There were 6 mice in the Responders group (animal numbers 1, 2, 4, 5, 6 and 8) and 4 mice in the Non-responders group (animal numbers 3, 7, 9 and 10) (Supplementary Fig. S4).

LC-MS/MS analysis of Verteporfin in mice blood samples

Verteporfin present in plasma samples of mice at necropsy was analyzed at the Proteomics facility of Cornell University (http://www.biotech.cornell.edu/brc/proteomics-and-mass-spectrometry). The LC-MS method was developed for detecting Verteporfin with high selectivity and sensitivity. We used targeted IDA method for identification of the analyte. There were two forms of Verteporfin (depending on the position of CH3 and H) found in standard as well as in the samples having the retention time 9.86 and 10.2 min respectively. Verteporfin was analyzed in control, Responder and Non-responder samples (n = 3 each).

Results

Transcriptome analysis of control and Verteporfin treated EMCA cells

We treated Type 1 EMCA cells (HEC-1-A and HEC-1-B) with VP at 10 nM for 3 h. This short time frame (half-life of VP is approximately 5–6 hours) and the lower concentration of VP was chosen to allow us to examine immediate responses to VP and to minimize confounding effects associated with severe damage and cell death. Under these conditions, VP causes several changes in gene expression within 3 h of treatment. We used RNA-seq to measure the transcriptomes altered by VP treatment. We generated 220–270 million reads per lane, filtered reads to have high quality scores and mapped 75–80% of those reads to the human genome. In order to get an overview over similarities and dissimilarities between samples, a hierarchical clustering to the heatmap function based on the sample distances was used (Figs 1 and S1). The strict counts of genes after VP-treatment tended to be less than in their non-treated counterparts in both HEC-1-A and HEC-1-B (Supplementary Fig. S2A). The differential expression of genes between the two groups was compared using a negative binomial test (Fig. 2). When considering the effect of VP treatment in both cell lines, 841 genes were upregulated and 251 genes down regulated after VP treatment (Fig. 1B). Treatment with VP was found to modulate many different genes and the top 20 upregulated (Table 1) and top 20 downregulated genes (Table 2) in EMCA cells after VP treatment were tabulated. The GO annotations of the differentially expressed genes for the combined cell line data were provided (Table 3). Positive regulation of TGF-β1 production, regulation of lipoprotein metabolic process, cell adhesion, endodermal cell differentiation, formation and development, and integrin mediated signaling pathway were among the significantly associated terms.

Table 1 Top 20 upregulated genes in HEC-1-A and HEC-1-B cells after VP treatment, ranked by logFC.
Table 2 Top 20 downregulated genes in HEC-1-A and HEC-1-B cells after VP treatment.
Table 3 Gene Ontology for top 20 upregulated genes in HEC-1A and HEC-1B after VP treatment (GO term results with FDR-adjusted p < 0.05), sorted by fold enrichment. The p-value was adjusted for False Discovery Rate (FDR).

When comparing DESeq-derived differentially-expressed genes by cell line, HEC-1-B showed the most marked change in RNA expression after VP treatment (Supplementary Fig. S2B, S2C). While many genes with low mean normalized count showed a high logFC after VP treatment, more genes with high mean normalized count showed general downregulation in HEC-1-B (Fig. 3 and Supplementary Tables S1, S2). For the analysis using only HEC-1-B samples, 8582 genes had differential expression with P < 0.05 after false discovery rate (FDR) correction. Of these, 3235 were downregulated and 5271 were upregulated in VP treated cells as compared to control cells. In HEC-1-B, differentially expressed genes with lowest p-values were VCAN, COL4A1, LAMC1, HSPA6, COL4A2, COL12A1, ITGAV, SPON1, KIF20A, ASPM. Notable was the fact that versican (VCAN), a gene that codes for a large extracellular matrix proteoglycan, had the lowest adjusted p-value (2.45e-216) and large absolute log2FC (log2FC = −3.463, lfcSE = 0.109) amongst all differentially expressed genes in HEC-1-B after VP-treatment. The tumor microenvironment has been found to contribute to VCAN mRNA expression (mostly in cleaved forms) in multiple tumor cell lines26,27.

Figure 3
figure 3

STRING network of protein-protein interactions between top 50 down-regulated and up-regulated genes in HEC1-B by adjusted p-value.

Gene co-expression

The VP-treated HEC-1-A cell line is more similar to VP-treated HEC-1-B in terms of transcriptomic gene expression; in fact, they are closer to each other than the control HEC-1-A and the control HEC-1-B (Figs 1 and 2). Functional enrichment analysis of differentially expressed genes in Type 1 cells after VP revealed that the extracellular matrix (ECM) organization Gene Ontology was the most significant (FDR p-value = 2.88e-15) (Table 3).

Protein product co-expression network

The resulting protein products of differentially expressed transcripts after VP-treatment were investigated for possible interaction with each other using STRING, which uses prior biological data (Figs 3 and S3, Supplementary Table S3)28. As the HEC-1-B cell line exhibited the most differential expression after VP-treatment, it was posited that the most differentially expressed genes would have strong interactions with each other. CDC23 and BUB1B, two genes crucially involved in mitotic checkpoint progression, were found to be the pair that had the best association from STRING among genes that were significantly differentially expressed in HEC-1-B after VP treatment (Supplementary Table S3). We obtained the Pfam protein domains and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways mostly associated with the differentially expressed genes using STRING (Fig. 3, Supplementary Table S5). Significantly enriched KEGG29 pathway protein products (Top 25) from genes differentially expressed in HEC-1-B after VP treatment were tabulated in Supplementary Table S4. Out of 124 genes in the KEGG cell cycle pathway (hsa-4110), 92 of them exhibited altered expression (adjusted p-value < 0.05) in HEC-1-B after treatment in VP; the top 25 are shown in Table 4. Enrichment in cell cycle genes was not present to the same degree in HEC-1-A. Further, we validated the RNASEQ data using qRTPCR in HEC-1-B cells treated with VP and observed that HSPA6 and C2orf88 genes were upregulated whereas CYR61, THBS1and ANKRD1 were downregulated after VP treatment. Based on the protein product co-expression data, we checked the expression of BUB1B, CDC23, MAD2L1, BUB3 and CDC27 genes and observed that these were downregulated after VP treatment in HEC-1-B cells. (Supplementary Table S6.) Similarly, we checked the expression of cell cycle genes CCRK, CDK2, CCND1, CCNE1 and E2F1 in both VP-treated EMCA cells (Supplementary Table S7) and VP-treated mice tumors (Supplementary Table S8) and observed that VP effectively inhibited the cell cycle genes. The downregulation of cell cycle proteins was established in VP-treated EMCA cells (Supplementary Fig. S5A). Further, we tested the role of YAP in regulating cell cycle proteins. YAP-knockdown in HEC-1-B cells by shYAP treatment resulted in downregulation of cell cycle proteins CCRK, CDK2 and Cyclin D1 (Supplementary Fig. S5B). These results show that VP is effective in inhibiting cell cycle both at transcription and translation levels. Interestingly, we also observed that VP downregulates pluripotency marker Oct4 either in vitro and in vivo (Supplementary Fig. S5C and S5D).

Table 4 Genes altered with an FDR adjusted p-value < 0.05 that belong to KEGG cell cycle pathways in HEC-1-B.

Effect of Verteporfin in vivo in nude mice

Since we observed significant cytotoxicity and induction of apoptosis by VP in EMCA cells, we employed a preclinical xenograft mouse model of EMCA to investigate the in vivo efficacy of VP on HEC-1-B cells in NCr nude mice. HEC-1-B GFP cells were injected into mice by subcutaneous (SC) administration. IP administration in this model achieved tumor regression (Fig. 4A,B) and importantly did not hasten time to euthanasia or decrease in weight between control and treated animals (Fig. 4C). These results show that VP is not toxic to mice. Based on the tumor volume data, we classified VP-treated mice into Responders and Non-responders (See Materials and Methods). Mean tumor volume in control mice was 558.34 mm3. Responders recorded mean tumor volume of 369.13 mm3 (33.89% decrease compared to control) whereas Non-responders had tumor volume of 619.15 mm3 (10.89% increase compared to control) (Supplementary Fig. S4A, S4C). Similarly, mean tumor weight of control mice was 0.302 g. In case of Responders, mean tumor weight decreased by 45.36% (0.165 g) compared to an increase of 18.54% (0.358 g) in Non-responders (Supplementary Fig. S4B). In order to check the concentration of VP in mouse blood, we analyzed their plasma samples using the LC-MS/MS method. Verteporfin was not detected in control samples and it was detected in Responder and Non-responder samples with good signal intensity. Responders had high amounts of VP (159.25 µg/ml) compared to non-responders (99.75 µg/ml) (Fig. 4D). These results indicate the bio-availability of VP and its effect on tumor regression. For further analyses, we used tumor samples of responders. Western analysis of control and VP treated tissues shows that VP is probably inducing tumor regression by acting upon cell cycle proteins (Fig. 4E). In accordance with our in vitro results, we observed a moderate effect of VP on EMCA tumor growth.

Figure 4
figure 4

(A) NCr nude mice showing tumor regression after VP treatments in a subcutaneous model after 21 days. n = 10/group. (B) Tumor volume growth curve and (C) Body weight curves in control and VP treated mice. (D) Quantitation of VP in mice plasma samples by LC-MS/MS analysis. (E) Western blot showing VP-induced inhibition of cell cycle proteins in subcutaneous tumors of mice. Error bars indicate Mean ± SEM.

Discussion

Verteporfin is a porphyrinic photosensitizer clinically used for the photodynamic treatment of age-related macular degeneration. YAP is a transcriptional co-activator and a potent oncogene30,31,32. Inhibition of the Hippo pathway leads to YAP activation, nuclear localization and cell proliferation in most cell types. The interaction of YAP with the TEAD family of transcriptional activators leads to DNA binding and transcription. YAP is amplified in many human cancers including breast, esophageal, hepatocellular, malignant mesothelioma, medulloblastoma and ovarian30,33,34,35,36. Based on the recent reports, VP is found to be a multi-target drug interacting with several proteins implicated in major cellular processes37. However, VP has been shown to inhibit autophagy in the absence of light activation both in vitro and in vivo38,39,40. Given the highly reactive nature of VP with light, Konstantinou et al.41 investigated the proposed mechanisms of non-light activated VP effects and suggested that VP-induced high molecular weight protein complexes (HMWC) require the presence of light. They observed that both singlet oxygen and radical generation mediate the formation of cross-linked oligomers and HMWC by VP in the presence of light. However, non-light activated, and light-independent VP effects have been demonstrated by several authors42,43,44. Verteporfin was identified as a potent inhibitor of cell growth in retinoblastoma cells, disrupting YAP-TEAD signaling and pluripotential marker Oct442. Verteporfin can also reverse the paclitaxel resistance induced by YAP over-expression in HCT-8/T cells without photoactivation through inhibiting YAP expression43. Cytotoxic effects of VP were mediated by light-independent production of reactive oxygen species and its anti-leukemic effects were also exerted in vivo without photoirradiation44. Similarly, Donohue et al. proposed that exposure to bright light could be harmful for VP-treated animals39. Based on these reports, in order to study the non-light activated effects of VP, we treated EMCA cells with VP for 3 h in the dark at 37 °C, followed by lysis with overhead fluorescence ambient laboratory lighting. Similar to previous findings, we observed maximum effect of VP in EMCA cells in inhibiting cell cycle and inducing apoptosis17. Based on the present results of RNA-seq analysis, VP induced modification of several genes belonging to different pathways, altering the transcriptomic etiology of EMCA cells.

Currently, VP has been widely used to understand the role of YAP1 in various cancers and is used for the treatment of various cancers17,45. Very few studies focused on the RNASeq analysis of EMCA pathophysiology. Mi et al. identified that the RACGAP1-STAT3-survivin signaling pathway is required for the invasive phenotype of uterine carcinosarcoma by applying RNA-Seq analysis to prospectively collected uterine carcinosarcoma tumor samples from patients46. Chen et al. studied the lncRNA transcriptome of endometrial cancers and adjacent normal endometrium from the same patients and compared with transcriptomes of other gynecologic malignancies including ovarian and cervical cancers47. The RNA-seq data of uterine corpus endometrial carcinoma (UCEC) samples identified potential genes associated with the development of UCEC. This group downloaded UCEC RNA-seq data from The Cancer Genome Atlas database48 and identified multiple key genes in UCEC and clinically relevant small molecule agents, thereby improving the understanding of UCEC and expanding perspectives on targeted therapy for this type of cancer. So far, there have been no reports on transcriptome of EMCA cells under in vitro conditions. In our study, we used RNA-seq and bioinformatic methods to better understand the EMCA cellular response to VP at the transcriptomic level, and possibly gain insight into the molecular mechanisms that might account for the pathophysiology of EMCA. Moreover, the underlying transcriptional response to VP suggests that EMCA cells carry out their function through different gene regulatory mechanisms. The results of the present study indicated that multiple differentially expressed genes are associated with EMCA pathophysiology. Differential expression analysis revealed changes in the expression patterns of many protein-coding genes previously reported to be involved in other cancers (Tables 34). Genes modulated by VP treatment include the genes related to positive regulation TGFβ1, lipoprotein metabolic processes, cell adhesion mediated by integrins, endodermal cell differentiation and extracellular matrix organization.

Even though we performed RNA-Seq analysis for both HEC-1-A and HEC-1-B cell lines, combining data related to both cell lines led to a result with not many significant pathways or gene-gene interactions. As in the PCA plot (Supplementary Fig. S1), the transcriptomic distance between HEC-1-A and HEC-1-B is too great to come up with many combined pathways. Hence, we focused our further analysis on HEC-1-B cell line data. Previously, we have shown that VP inhibits cell cycle progression of Type 1 EMCA cells17, hence we analyzed KEGG cell cycle pathways in HEC-1-B. Out of 124 genes in the KEGG cell cycle pathway, 92 of them exhibited altered expression after VP treatment. These results also corroborate with our in vivo studies, as we observed that tumor regression in mice is induced by inhibition of cell cycle proteins. Amongst all genes that were differentially expressed in HEC-1-B after VP treatment, naïve of biological annotation, mRNA transcripts of VCAN were most likely to be downregulated with VP treatment (Supplementary Fig. S2). The role of Versican in cell adhesion, migration, and proliferation has been extensively studied49,50, and in cancer, stromal expression of VCAN was strongly correlated with disease recurrence26. The downregulating effect of VP treatment on VCAN expression for the HEC-1-B cell line is promising, but the precise mechanism of this effect remains to be elucidated, especially considering the fact that a contrary upregulating effect on VCAN expression was found in HEC-1-A, or when both HEC-1-A and HEC-1-B were analyzed together (Table 1). Cyclin-dependent kinase 20 (CDK20) or CCRK, is a member of CDK family with strong linkage to human cancers. Recent studies reported the consistent overexpression of CCRK in cancers arising from brain, colon, liver, lung and ovary. The signaling molecules perturbed by CCRK are divergent and cancer-specific, including the cell cycle regulators CDK2, cyclin D1, cyclin E and RB in glioblastoma, ovarian carcinoma and colorectal cancer and lung cancer51. Overexpression of CCRK increases cyclin D1 expression, which suggests the role of CCRK in the control or cell proliferation via regulation of cyclin D1 expression52. Recently Zapiecki et al. showed that significantly higher expression of cyclin D1 and cyclin E was detected in patients dying from endometrial cancer53. Cyclin D1 plays an important role in LSD1-regulated estrogen-induced endometrial cancer cell proliferation54. This group also verified the positive correlation between LSD1 and cyclin D1 in endometrial cancer tissues by immunohistochemistry. In conformity with previous findings, our studies also indicate that inhibition of cyclin D1, CCRK and CDK2 are involved in inhibition of growth and proliferation in EMCA cells and tumor regression in mice. Normal adult stem cells and cancer stem cells maintain expression of Oct3/4, consistent with the stem cell hypothesis of carcinogenesis. Hence, a strategy to target “cancer stem cells” is to suppress the Oct4 gene expression55. Recently Ding et al. studied the characteristics of CD133+ cells isolated from endometrial cancer56. They proposed that CD133+ cells increased expression of embryonic stem cells markers including Oct4, Nanog, Sox2 than CD133 cells. Verteporfin was shown to downregulate Oct4 expression in retinoblastoma cells41. Our findings corroborate with previous results under both in vitro and in vivo conditions and also suggest a mechanism of action of VP.

The administration of VP in subcutaneous mouse model identified both Responders and Non-responders to this treatment. Our results show that VP is not toxic to mice without significant effect on their body weight. Visudyne®, the commercial liposomal formulation of VP was developed for photodynamic therapy. Based on the calculations of low verteporfin-to-lipid ratio of Visudyne, Donohue et al. administered 20 mg/kg VP and observed that VP tumor accumulation following i.v. administration was very poor in pancreatic ductal adenocarcinoma (PDAC) mouse model39. To increase VP tumor accumulation, they developed an alternate micellar formulation with DSPE-mPEG2000 dissolved in PBS. They observed that VP inhibited autophagy in vivo but did not reduce tumor volume or increase survival as a single agent in PDAC mouse model. They also observed that VP in combination with gemcitabine moderately reduced tumor growth and enhanced survival compared to gemcitabine alone. Similar to their findings Zhao et al. demonstrated that VP in combination with a pan-RAF inhibitor LY3009120 significantly enhanced tumor regression in KRAS-mutant pancreatic cancer57. Our results support their findings as we observed that VP failed to show antitumor activity in non-responders, even though the drug accumulated at the tumor site based on our LC-MS/MS analysis. These results warrant further standardization of solubility of VP, dosing schedule of VP, route of administration, as well as use of VP in combination with other chemotherapeutic agents for effective treatment in EMCA. On-going experiments in our laboratory include development of orthotopic model of EMCA, administration of VP as mPEG-PLGA nanoparticles and IV administration of VP with pharmacokinetic studies. The present study revealed multiple key genes of pathological significance in EMCA as well as mechanism of activity of VP with the ultimate goal of improving our understanding of molecular profile of EMCA which expands our perspective on targeted therapy for this cancer. Taken together our results suggest that VP may be repurposed for the treatment of advanced endometrial cancer.