PC-3 (cat. no. CRL-1435), LNCaP (cat. no. CRL-1740), DU145 (cat. no. HTB-81) and 22Rv1 (cat. no. CRL-2505) prostate cancer cell lines were purchased from the American Type Culture Collection. Cells were cultured in RPMI 1640 (PC-3, LNCaP and 22Rv1) or DMEM (DU145) containing 10% fetal bovine serum, penicillin (100 U/ml) and streptomycin (100 µg/ml). All cell culture reagents were obtained from HyClone (Cytiva). All other reagents not specified were supplied by MilliporeSigma. All cell lines tested negative for mycoplasma contamination using the Mycoplasma PCR Detection kit (Intron Biotechnology, Inc.). These cell lines have been authenticated in the three years using short tandem repeat analysis.

Male Balb/C nude mice (4–5 weeks old; 18–20 g; n=10) were purchased from Charles River Laboratories Japan, Inc. Mice were housed in laminar flow cabinets under specific pathogen-free condition (37°C; 12-h light/dark cycle; 60% relative humidity and free access to food and water). Mice were anesthetized by the inhalation of isoflurane (Terrell; Piramal Critical Care Inc.) in oxygen (2 l/min): Induction with 4% isoflurane for 2 min in an anesthetic induction chamber, followed by maintenance with 2% isoflurane for 1–2 min after transferring to a nose cone. The 22Rv1 cells (5×10⁶ cells; 100 µl of cell suspension) were subcutaneously injected into the right flank of each mouse. The humane endpoints were when the largest tumor size was >20 mm in diameter. None of the mice reached the endpoints of the present study. When the tumor reached ~180 mm³, mice were randomly divided into two groups (five in each group) and intraperitoneally injected with 20 mg/kg CsA or DMSO, every other day for 14 days. After 14 days, mice were humanely sacrificed under overdosed isoflurane. Mice were placed into a chamber filled with vapor of the anesthetic isoflurane until respiration ceased (within 2 min) and the tumors were excised. The tumor volume was calculated using the following formula: V=(L × W²) × 0.5 (V, volume of tumor; L, length of tumor; W, width of tumor). All animal experiments were performed in accordance with protocols approved by the Institutional Animal Care and Use Committee of the Asan Institute for Life Sciences at the ASAN Medical Center, University of Ulsan College of Medicine, Seoul (approval no. 2021-13-234).

After cells were transfected with the short interfering (si)RNAs for 48 h, MTT assay was performed to assess cell growth according to the manufacturer's instructions (MilliporeSigma). The assay was quantitated by measuring the absorbance at 570 nm on a BioTek SynergyMx microplate reader (BioTek Instruments, Inc.).

PC-3 (2–5×10³ cells/well), LNCaP (8×10³ cells/well), DU145 (1–3×10³ cells/well) and 22Rv1 (2–5×10³ cells/well) cells were seeded into 6-well plates. Cells were treated with CsA (10 µM) or siRNAs for 9 days (once every 3 days). The cells were fixed with 3.7% paraformaldehyde and stained with 0.5% crystal violet (MilliporeSigma) for 15 min at room temperature. The number of colonies, defined as >50 cells/colony, was counted using ImageJ (1.8.0 172; National Institutes of Health).

PC-3 cells were treated with 10 µM CsA in 0.1% ethanol (vehicle) for 24 h. Microarray experiments were performed as we previously described (18,19). The microarray data are available through the Gene Expression Omnibus database (GSE109505; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE109505). Rigorous data preprocessing and single channel array normalization (SCAN) were performed and microarray probes were mapped to gene symbols as previously described (20–22). Of the 20,661 mapped genes, 17,629 protein-coding genes were selected for further analyses. The internal clusters were validated through hierarchical clustering and principal component analysis (PCA) (20,21).

The GSE67157 and GSE109505 datasets were used to construct the PC-3 cell-specific interactome. Rigorous preprocessing, such as quality control testing, normalization and batch effect adjustment, was conducted as previously described (21,23,24). Of the common 11,506 genes from two datasets, 5,000 genes of high variance were selected for implementing the Algorithm for the Reconstruction of Gene Regulatory Networks (ARACNe) (20,21). The list of human transcription factors (TFs) was obtained from the Animal Transcription Factor Database 2.0 (AnimalTFDB 2.0) (23) and used for Master Regulator Analysis (MRA). ARACNe preprocessing and MRA-Fisher's exact test (FET) analysis were run in geWorkbench software version 2.6.0 (http://wiki.c2b2.columbia.edu/workbench/i-ndex.php/Home) as described in our previous reports (20,21,25,26). GSE3325 and GSE35988 data were used to assess the changes in gene expression profiles during prostate cancer progression. A more detailed description is provided in Fig. S1 and its legend.

SAM was used to identify differentially expressed genes (DEGs) from the GSE67157, GSE109505, GSE3325 and GSE35988 data (27). A tuning parameter, delta of 0.4, was optimized to give the cutoff for significance with the estimation of the false discovery rate (FDR) q-value threshold of 0.01. GSEA (Hallmark Gene Set from the Molecular Signature Database) was performed to obtain a biological interpretation of clinically significant CsA-specific DEGs (28). Universal concept signature scores were calculated for E2F8 target genes and concept signature enrichment analysis performed for deep functional assessment of the pathways. A more detailed description is given in Fig. S1 and its legend.

PC-3, LNCaP, DU145 and 22Rv1 cells were treated with CsA at the indicated concentrations for various times. Total RNA was extracted using the RNeasy Mini Kit (Qiagen GmbH) and reverse transcribed with SuperScriptIII First-Strand (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturers' instructions. RT-PCR and real-time PCR were performed using specific primers for E2F8 (16,21), MELK (29), β-actin (30) or 18S (31). The primer sequences are provided in Table SI. The thermocycling conditions for PCR were: 95°C for 5 min, followed by 31 cycles of 95°C for 45 sec, 51°C (MELK and β-actin) or 54°C (E2F8) for 45 sec, and 72°C for 45 sec. PCR products were separated by 1% agarose gel electrophoresis and visualized with SYBR Safe DNA Gel Stain (Invitrogen; Thermo Fisher Scientific, Inc.). Agarose gel electrophoresis images were acquired with Gel Documentation XR System (Bio-Rad Laboratories, Inc.). Real-time PCR was conducted using StepOne Real-Time PCR system (Thermo Fisher Scientific Inc.).

For western blotting, the crude extracts were prepared by incubation with RIPA buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1% triton X-100, 0.5% sodium deoxycholate, 0.1% SDS, 0.5 M EDTA) containing protease and phosphatase inhibitor cocktails (MilliporeSigma). The protein concentrations were determined by BCA assay kit (Thermo Fisher Scientific, Inc.). The samples (10–30 µg for each) were resolved in 7% (E2F8), 8% (MELK), or 10% (E2F1 and β-tubulin) SDS-PAGE gels and transferred onto NC membranes (Bio-Rad Laboratories, Inc.). The membranes were blocked using 5% w/v skimmed milk (BD Difco; BD Biosciences) in Tris-buffered saline containing 0.1% Tween-20 (TBST) for 1 h at room temperature. The membranes were probed with the anti-E2F8 (1:2,000; cat. no. A303-039A; Bethyl Laboratories, Inc.), anti-MLEK (1:1,000; cat. no. 2274S; Cell Signaling Technology, Inc.), anti-E2F1 (1:500; cat. no. sc-193; Santa Cruz, Inc.), and anti-β-tubulin (1:5,000; cat. no. T4026; MilliporeSigma, Inc.) antibodies for 1 h at room temperature. After washing three times with TBST, the membranes were incubated with a goat anti-Rabbit IgG-HRP antibody (1:5,000; cat. no. A120-101P, Bethyl Laboratories, Inc.) or a goat anti-Mouse IgG-HRP antibody (1:5,000; cat. no. A90-116P, Bethyl Laboratories, Inc.) for 1 h at room temperature. The signals were determined by the enhanced chemiluminescence reaction (ECL; Amersham; Cytiva). X-ray films were scanned and analyzed using EPSON Scan Software (EPSON Expression 11000XL, Seiko Epson Corporation). The data shown are representative of at least four independent experiments. Full scan images are shown in Appendix S1.

Cells were transfected with 50 nM siControl (21) or siRNAs against E2F8 [siE2F8-1 and siE2F8-2 (21)] or MELK [siMELK-1 (29) and siMELK-2 (32)] for 48 h at 37°C using Lipofectamine^® RNAiMAX reagent (Invitrogen; Thermo Fisher Scientific, Inc.). The siRNAs were obtained from Genolution Pharmaceuticals Inc. (Seoul, Korea). The siRNA sequences are provided in Table SII. The time interval between transfection and subsequent experiments was 48 h or 72 h for MTT assay, 9 days for colony formation assay (see also Figure legends).

The Cancer Genome Atlas (TCGA)-Prostate Adenocarcinoma (PRAD) gene expression data were obtained from Xena Functional Genomics Explorer (https://tcga.xenahubs.net/download/TCGA.PRAD.sampleMap/HiSeqV2.gz). The E2F8 gene signature, which was inferred in MRA-FET analysis, was extracted and the continuous value of expression levels in each gene was discretized to equal-width bins (8 bins) using the unsupervised discretization method available in the Information-Theoretic Measures (Infotheo) R version 4.1.1 (http://www.R-project.org/) (11). The fast greedy equivalence search (FGES)-discrete algorithm was used to identify E2F8-interacting genes (33). A more detailed description is given in Fig. S1 and its legend.

The gene expression data of TCGA-primary tumor samples were downloaded using the TCGA Biolinks package (ver. 2.14.1, http://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html). Before normalization, duplicated FFPE samples were removed. Gene reads were normalized in counts per million (CPM) using the edgeR package (ver. 3.28.1, http://bioconductor.org/packages/release/bioc/html/edgeR.html) and log₂ transformation (CPM+1). The Ensembl gene IDs were mapped into HGNC symbols using the biomaRt package (ver. 2.42.1, http://bioconductor.org/packages/release/bioc/html/biomaRt.html) and the expression data of androgen receptor (AR), MELK and E2F8 extracted for coexpression analysis.

TCGA-PRAD RNA-seq raw bam files (n=554) were downloaded from NCI Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/). Patient response data for chemotherapy or radiotherapy were provided by Xena Functional Genomics Explorer (https://xenabrowser.net/hub/) GDC Hub (34). When the therapeutic response data were mapped to the available PRAD RNA-seq samples, 313 patients were identified as having a complete response (CR), 34 patients as partial response (PR), 29 patients as progressive disease (PD) and 27 patients as stable disease (SD). Integrated System for Motif Activity Response Analysis (ISMARA, http://ismara.unibas.ch/mara/) was used to infer regulatory networks and sample-specific TF motif activities from gene expression data (35). ISMARA allowed the mapping of transcriptomic profiles to a lower-dimensional inferred TF activity space, largely preserving the relationships between samples. The TF activities driving expression changes were calculated and the TF activity differences between therapeutic response and non-response patients were examined. The TF activities significantly associated with therapy response were defined by P-value <0.05 and absolute mean activity difference >0.008.

The Kaplan-Meier survival curve and log-rank test were used to determine overall survival curves as previously described (36,37). Using median gene expression values from cBioPortal transcriptomic data as a bifurcating point, the samples were divided into high- and low-expression groups and the survival rates were compared between the two groups. The Cox proportional hazards model was applied to estimate hazard ratios (HRs) and 95% confidence intervals (CIs). Pearson correlation was used to determine the correlations of the expression levels among E2F8, MELK, or AR from the TCGA-PRAD and cBioPortal data. A comparison of means among experimental groups was performed using one-way ANOVA followed by Bonferroni's multiple comparison test. P<0.05 was considered to indicate a statistically significant difference.

To determine the antitumor activity of CsA in prostate cancer, MTT and colony formation assays were first performed. CsA suppressed cell growth and colony formation in PC-3, LNCaP, DU145 and 22Rv1 prostate cancer cells (Fig. 1A and B). In addition, CsA inhibited tumor volume and weight in 22Rv1 cell xenograft mouse model (Fig. 1C and D). These results demonstrate that CsA has an anticancer activity against prostate cancer.

To understand the antitumor mechanism of action of CsA, microarray experiments we performed using CsA-treated PC-3 cells (GSE109505). Hierarchical clustering analysis and PCA showed that CsA-treated and untreated cells were clustered into two discrete groups (Figs. S1, S2A and S2B), indicating that CsA induces a distinct change in gene expression profiles. Analysis using SAM found that CsA significantly affected the expression levels of 3,319 genes (Fig. S2C). Among these 3,319 differentially expressed genes (DEGs), 2,500 (593+1,907) genes were downregulated in CsA-treated cells, whereas 819 (278+541) genes were upregulated in CsA-treated PC-3 cells (Fig. S2D).

To assess the clinical significance of the 3,319 DEGs, they were compared with the DEGs derived from the transcriptomic data of patients with prostate cancer (GSE3325). First 3,654 DEGs between metastatic and nonmetastatic (benign or primary) cancer were obtained: 1,398 genes were upregulated and 2,256 genes were downregulated in patients with metastatic cancer (Fig. S2D; leftmost circles). It was identified that 871 (593+278) DEGs showed anti-similarity (or inverse correlation) between the 3,319 DEGs from GSE109505 and the 3,654 DEGs from GSE3325 (Fig. S2C and D). These 871 anti-similar DEGs were termed the ‘clinically significant CsA-induced gene expression (CCI) signature’. GSEA revealed that this CCI signature is associated with 19 hallmark pathways (Fig. S2E). In particular, CsA activates cell death-related pathways and inhibits cell cycle-related pathways, providing the mechanistic explanation for the anticancer activity of CsA (Fig. 1).

To elucidate the molecular mechanism by which CsA regulates the CCI signature, the prostate cancer cell-specific transcriptional interactome was analyzed using ARACNe and MRA algorithms. These identified 28 transcription factors as master regulators (MRs) that control the CCI signature (Fig. 2A and Table I): 27 MRs were markedly downregulated in CsA-treated PC-3 cells, whereas only one MR (PHF1) was upregulated (Fig. 2A; heatmap). Overall, the expression patterns of these MRs correlated with those of their target genes (Fig. 2A; red bars). Based on ‘markers in the intersection set’ (Table I), E2F8 was identified as a top MR that affected the majority of genes belonging to the CCI signature (216 target genes out of 871; Figs. 2A and S3). These 216 target genes of E2F8 were mainly associated with cell cycle or proliferation pathways (Fig. S3). CsA markedly reduced E2F8 expression (Fig. 2A), which correlated positively with the expression levels of the most target genes (bar plot in Fig. 2A). These results suggested that E2F8 mainly acts as a transcriptional activator.

To assess the clinical relevance of the MRs, the data of patients with prostate cancer (GSE3325 and GSE35988) were analyzed. Among these 28 MRs, the expression levels of 10 MRs were markedly upregulated during prostate cancer progression (E2F8 in Fig. 2B and other nine MRs in S4). Kaplan-Meier analysis showed that high expression levels of these 10 MRs are associated with worse prognosis in patients with prostate cancer (E2F8 in Fig. 2C and other nine MRs in S5). In particular, E2F8 showed the highest hazard ratio (HR=3.028; P=0.0002; Figs. 2C and S5B). These results indicated that E2F8 acts as a clinically significant MR crucial for governing a large portion of the CCI signature. In addition, the findings suggested that E2F8 serves as a useful prognosis marker for prostate cancer.

To confirm whether CsA suppressed E2F8 expression, RT-PCR, real-time PCR and western blot analysis we performed in various types of prostate cancer cells. CsA inhibited E2F8 mRNA (Figs. 3A and B and S6A) and protein expression (Fig. 3A and B) in a time- and concentration-dependent manner. As PC-3 and DU145 cells do not express functional AR, the results suggested that CsA suppressed E2F8 expression through an AR-independent mechanism.

To assess the role of E2F8 in prostate cancer, the effect of siRNAs against E2F8 (siE2F8-1 and −2) were examined using MTT and colony formation assays. The siE2F8s effectively inhibited E2F8 expression in all tested cell lines (Fig. S6B). It was found that E2F8 knockdown suppressed cell growth and colony formation in all tested cells (Fig. 3C and D). These results suggested that E2F8 represents an attractive therapeutic target for prostate cancer.

To elucidate the molecular mechanism by which CsA suppresses E2F8 expression, Bayesian network analysis was employed, which is an effective approach to model causal relationships between observed biological data and gene expression data (38,39). The FGES-discrete algorithm identified MELK as an E2F8-interacting gene. Analysis of TCGA-PRAD and cBioPortal transcriptomic data showed the strong positive correlation between E2F8 and MELK expression in patients with prostate cancer (r=0.7704, P<0.0001 and r=0.7415, P<0.0001, respectively; Fig. 4A). In addition, analysis of the GSE3325 data showed that MELK is markedly upregulated during prostate cancer progression (Fig. 4B). Whether MELK causally regulates E2F8 expression was then examined. The siRNA-based knockdown of MELK (siMELK-1 and −2) markedly reduced E2F8 expression in all tested cells (Fig. 4C). In contrast to E2F8, siMELKs did not affect the expression of E2F1, one of the identified MRs (Fig. 2 and Table I), which plays a crucial role in prostate cancer growth (20,40). Altogether, these results indicate that MELK is a crucial upstream signaling molecule for controlling E2F8 expression in prostate cancer cells.

To determine whether CsA suppresses MELK expression, RT-PCR, real-time PCR and western blot analysis were performed in various prostate cancer cells. CsA inhibited MELK expression at the protein level (Fig. S7A), but not at the mRNA level (Fig. S7A and B), suggesting that CsA regulated MELK at the posttranscriptional or translational level. To assess the functional importance of MELK in prostate cancer cells, whether forced suppression of MELK expression affects prostate cancer growth was investigated using MTT and colony formation assays. The siMELKs inhibited cell growth and colony formation in various prostate cancer cells (Figs. 4D, S7C and S7D). Altogether, these results indicate that the MELK-E2F8 signaling axis plays a crucial role in prostate cancer biology.

AR is a crucial driver of CRPC progression and its expression is frequently upregulated during prostate cancer progression (41–44). The present study showed that E2F8 expression increases during prostate cancer progression and exerts oncogenic activity (Figs. 2 and 3). Based on these observations, the interplay between AR and E2F8 signaling was examined. The analysis of TCGA-PRAD and cBioPortal transcriptomic data found that there was no significant coexpression correlation between E2F8 and AR or between MELK and AR expression levels (Fig. 5A and C).

Nonetheless, Kaplan-Meier analysis showed that overall survival was greatly reduced in patients with prostate cancer with high levels of both E2F8 and AR expression, compared with those with low levels of both genes (HR=4.33, P=0.0019; Fig. 5B). Similar results were obtained from patients with prostate cancer with high levels of both MELK and AR expression, compared with those with low levels of both genes (HR=3.3, P=0.0041; Fig. 5D). Altogether, these results suggest that the MELK-E2F8 axis and AR signaling act independently, but their additive effects worsen the clinical outcomes of prostate cancer. Hence. simultaneous inhibition of E2F8 and AR could be a promising therapeutic strategy for patients with prostate cancer with concurrent overexpression of these molecules.

To investigate the translatability of the findings, LNCaP cells (AR-positive, AR-V7-negative cells) were treated with siE2F8 and the AR antagonist MDV3100. The treatment with siE2F8 increased the sensitivity of prostate cancer cells to MDV3100 (Fig. 5E), suggesting a potential therapeutic strategy for improving the efficacy anti-androgen therapy. In addition, siE2F8 had a similar effect as siMELK in enhancing the efficacy of MDV3100 (Fig. 5F). Finally, the ISMARA algorithm showed that E2F8 activity is markedly higher in patients with SD than in those with CR (Fig. 5G), suggesting that E2F8 is a predictive marker for therapeutic response in patients with prostate cancer.