Association of OIP5 upregulation with pRCC tumorigenesis and progression
OIP5 was reported to be a component gene in a multigene set predicting the risk of prostate cancer recurrence (38); its upregulation associates with adverse features in ccRCC and bladder cancer (19, 29), supporting a general involvement of OIP5 in urogenital cancers. To investigate this possibility, we examined OIP5 expression in pRCC using a tissue microarray (TMA) containing 40 pairs of pRCC and 74 pairs of ccRCC tumors with the adjacent non-tumor kidney (AJK) tissues from 20 and 37 patients, respectively. In comparison to AJK tissues, pRCC tumor tissues expressed a significant OIP5 upregulation (Fig. 1a, b). OIP5 expression was further increased in advanced T stage tumors (Fig. 1b). Consistent with a previous report, OIP5 upregulation occurred in Grade 2-3 ccRCC tumors compared to the AJK tissues; nonetheless, we could not demonstrate OIP5 upregulation in Grade 1 ccRCC compared to the AJK tissues (Additional file 1: Fig. S1), suggesting a role of OIP5 in ccRCC progression. By using the TCGA RNA-sequencing data organized by UALCAN (ualcan.path.uab.edu/home) (37), OIP5 upregulation at the mRNA level in pRCC tissues was observed (Fig. 1c); the upregulations reflects the level of severity and the order of unfavorable outcome of pRCC with higher expression levels in T2P over T1P tumors, CIMP tumors over other subtypes (Fig. 1d), stage 3 tumors over stages 1-2 tumors, stage 4 over stage 3 tumors (Fig. 1e), and N1 (lymph node metastasis) over N0 tumors (Fig. 1f). Consistent with its associations with adverse tumor features, OIP5 expression robustly stratifies pRCC tumors into a high- and low-risk group based on overall survival possibility (Fig. 1g). Among the 10 patients in the OIP5-high group, 7 died in a rapid time course (Fig. 1g). Collectively, these observations support a strong association of OIP5 with pRCC tumorigenesis and progression.
OIP5-mediated enhancement of pRCC tumorigenesis along with network alterations
Attributed to the uncommon status of pRCC, there are only limited number of confirmed pRCC cell lines available. ACHN is the most widely used and confirmed metastatic pRCC cell line; the cells have the typical feature of c-MET polymorphism detected in pRCC (39, 40). ACHN is likely the only confirmed metastatic pRCC cell line (39). To analyze the functional impact of OIP5 on pRCC tumorigenesis, we stably expressed OIP5 in ACHN cells (Fig. 2a). In comparison to ACHN EV (empty vector) cells, ACHN OIP5 cells displayed elevated abilities for proliferation (Fig. 2b), colony formation (Fig. 2c, Additional file 2: Fig. S2a), invasion (Fig. 2d, Additional file 2: Fig. S2b) and growth in soft agar (Fig. 2e, Additional file 2: Fig. S2c). We have also established the EV and OIP5 stable lines in the commonly used 786-O ccRCC cells, and OIP5 overexpression did not affect all of the above oncogenic events observed in ACHN cells in vitro (data not shown), which suggests a certain level of specificity of OIP5 in promoting pRCC. In vivo, OIP5 enhanced the growth of ACHN cell-produced xenografts compared to tumors produced by ACHN EV cells (Fig. 2f); mice bearing ACHN OIP5 tumors reached endpoint faster compared to animals with ACHN EV cell-produced tumors (Fig. 2g). The overexpression of OIP5 in ACHN OIP5 tumors was confirmed (Additional file 3: Fig. S3a). The ACHN OIP5 tumors show a significant increase of CDK2 expression largely in the nuclei (Additional file 3: Fig. S3b); the functions of this are not clear as no upregulations of the relevant cyclins (cyclin A and cyclin E) was observed (data not shown).
To further analyze factors and networks utilized by OIP5 in enhancing ACHN cell-produced xenografts, RNA-sequencing (RNA-seq) was performed on ACHN EV and ACHN OIP5 tumors at three per group. Gene set enrichment analysis (GSEA) was conducted on differentially expressed genes obtained in the setting of OIP5 vs EV. When enrichment in the oncogenic gene sets (C6) collection was analyzed using FGSEA (fast gene set enrichment analysis), we observed that genes downregulated (DN) in cells with activation (UP) of ERB2, MEK, and mTOR were also downregulated in ACHN OIP5 tumors compared to ACHN EV tumors (Fig. 3a), suggesting OIP5 suppressing those genes that are downregulated by ERB2, MEK, and mTOR. Similarly, ACHN OIP5 tumors also display downregulation of EGFR-downregulated genes (Additional file 11: Table S1). The serine/threonine kinase 33 (STK33) is a synthetic lethal interacting protein of KRAS mutant, i.e. cells expressing KRAS mutant rely on STK33 for survival (41). Knockdown of STK33 in acute myeloid leukemia cells led to upregulation of a set of genes (STK33-UP) (41), suggesting a potential inhibition of these genes by STK33. These gene expressions were also reduced in ACHN OIP5 tumors (Fig. 3a, Additional file 11: Table S1). To test the reliability of the enrichment obtained by FGSEA, GSEA was further conducted using a more stringent platform: EGSEA. Ensemble gene set enrichment analysis produces a consensus gene set ranking (enrichment) with the combination of multiple (up to n=12) algorithms (42). With the maximal stringent condition using all 12 algorithms, EGSEA revealed within the top 12 ranks the downregulation of the ERB2- and MEK-suppressed gene sets in ACHN OIP5 tumors (Additional file 4: Fig. S4); downregulation of genes in cells with STK33 knockdown was observed in multiple setting (Additional file 4: Fig. S4) which is consistent with the enrichments derived from using FGSEA (Additional file 11: Table S1). All top 12 ranked gene sets obtained by EGSA (Additional file 4: Fig. S4) are also included in those produced by FGSEA (Additional file 11: Table S1). It is intriguing that VEGFA-suppressed genes in HUVEC (human umbilical vein endothelial cell) cells were also downregulated in ACHN OIP5 tumors (Additional file 4: Fig. S4, Additional file 11: Table S1). Based on the overall gene set enrichment within the oncogenic gene set (C6, MSigDB) collection (Additional file 11: Table S1), we can summarize that in ACHN OIP5 xenografts, the RB pathway is inhibited and the signaling processes of STK33, BMI1, EZH2, MYC, WNT, VEGFA, and EGFR/ERB2 are enhanced (Fig. 3b).
We further examined gene set enrichment within the Hallmark gene set collection using FGSEA. The analyses revealed downregulations of inflammatory response, TNFα_via_NFκB signaling (NFκB-regulated genes in response to TNFα), and complement gene expression (Hallmark_Component, normalized enrichment score/NES: -1.48, padj 0.013) (Additional file 5: Fig. S5a; Additional file 12: Table S2). Additionally, ACHN OIP5 xenografts exhibited upregulations in gene sets regulating fatty acid metabolism and cholesterol homeostasis (Fig. 3C; Additional file 12: Table S2). These enrichments were also produced by EGSEA (Additional file 5: Fig. S5b, c). Several processes are enhanced in ACHN OIP5 tumors, which include oxidative phosphorylation, the expression of E2F and MYC targets, EMT (epithelial mesenchymal transition), mTORC1 signaling, and adipogenesis (Additional file 12: Table S2). Enrichment in glycolysis in ACHN OIP5 tumors was obtained by FGSEA (Additional file 12: Table S2), which was also confirmed by KEGG pathway analysis using EGSEA (Additional file 6: Fig. S6). Evidence thus suggests a metabolic switch to Warburg metabolism in ACHN OIP5 tumors.
Association of OIP5-related differentially expressed genes with CIMP subtype
In comparison to other pRCC subtypes, CIMP tumors have a Warburg metabolic shift (12), indicating an association between OIP5-affected genes and CIMP. This notion is supported by the elevation of OIP5 expression in CIPM pRCC tumors (Fig. 1d). To investigate this possibility, we firstly defined the differentially expressed genes (DEGs) in ACHN OIP5 tumors vs ACHN EV tumors as those with p.adj < 0.05 and fold change > |1.5|; a total of 1128 DEGs were derived (Additional file 13: Table S3). In these DEGs, the top upregulated genes include WNT7A, FGF1, CNTN1,(43) SOX2, and others, which are known for their facilitative roles in tumorigenesis. The top 20 clusters enriched in these DEGs contain those that regulate urogenital system development, blood vessel morphogenesis, hippo pathway, cell surface receptor signaling, pathway in cancer, epithelial cell proliferation, and others (Fig. 4a, Additional file 14: Table S4). Individual terms in these enriched clusters forms a network connection (Additional file 7: Fig. S7a). These pathways are clearly relevant to tumorigenesis. DEGs are clustered in ACHN OIP5 tumors vs ACHN EV tumors (Fig. 4b).
To confirm the relevance of these DEGs derived from ACHN cell-produced xenografts in pRCC pathogenesis, we analyzed their relationship to DEGs derived from primary pRCCs relative to OIP5 expression. In the TCGA Pancancer pRCC dataset within cBioPortal, high OIP5 expression robustly separates pRCC tumors into a high and low risk group based on their overall survival (OS) possibilities (Fig. 1g). From these two groups, we obtained 873 DEGs defined by q < 0.05 and fold change ≥ |2| (Additional file 15: Table S5). These primary pRCC-derived DEGs share 66 overlap DEGs (Overlap66) with the xenograft-derived DEGs (Table 1, Fig. 4c). The alterations in their expressions in normal kidney tissues (n=30) and pRCC tumors (n=290) at different stages are presented in Additional file 8: Fig. S8. The genes with further elevations in Stage 3-4 tumors include SLC7A11, PCSK5, STC2, PLK1, TK1, TRIB3, and SRXN1 (Additional file 8: Fig. S8).
Among these 66 DEGs, 8 and 41 genes are not known for associations with cancer and ccRCC respectively (Additional file 16: Table S6a); only PLK1 was reported to be a component gene in a prognostic multigene of pRCC (Additional file 16: Table S6a). Overlap66 is novel to pRCC. 46 out of 66 DEGs significantly predict overall survival (OS) possibility with some being individually efficacious based on their p values: 6.73e-10 for PCSK5, 1.3e-10 for C116orf75 (RMI2), 1.84e-13 for ATAD2, and others (Additional file 16: Table S6a). Furthermore, 33 DEGs retain their predictive significance after adjusting for age at diagnosis, sex, and T stage (Additional file 16: Table S6a).
The potentials of the 33 DEGs as prognostic biomarkers are in accordance with their expression status in CIMP. C116orf75, SRXN1, TK1, and TRIB3 positively predict poor OS (Table 1, Additional file 16: Table S6a); they are notably upregulated in CIMP tumors (Fig. 4d). Reversely, CCDC106, CX3CL1, LYNX1, and SPATA18 are negatively associated with poor OS (Table 1, Additional file 16: Table S6a); their expressions are particularly downregulated in CIMP tumors (Fig. 4e). In all 46 genes with their expression associated with OS shortening, 11 show no alterations in gene expression in CIMP tumors (Table 1); for the remaining 35 genes, their positive and negative predictions of OS shortening correlate with their respective upregulation and downregulation in CIMP tumors (Table 1). This correlation of expression was not observed in tumors vs non-tumor tissues (Table 1). In view of CIMP tumors having the poorest OS possibility,(12) the association of these gene expression with CIMP tumors supports their potential as prognostic biomarkers.
Robust prognostic biomarker potential of Overlap66 and its sub-multigene panels
Following the above analyses, we examined the OS-related prognostic potential of Overlap66 as a multigene panel. The expression data for these DEGs along with the relevant clinical data were retrieved from the Pancancer pRCC dataset within cBioPortal. Risk scores for individual tumors were calculated as ∑(coefi x Geneiexp)n (coefi: Cox coefficient of genei, Geneiexp: expression of Genei, n=66). Coefs were obtained using the multivariate Cox model. Overlap66 risk scores efficiently predict OS shortening using both univariate (UV) and multivariate (MV) Cox models (Additional file 7: Fig. S7b). The MV model consists of the risk scores, age at diagnosis, sex, and T stage (Additional file 7: Fig. S7b). With the cutoff point optimized using the Maximally Selected Rank Statistics (Additional file 9: Fig. S9), Overlap66 effectively stratifies the risk of fatality (possibility of OS) and relapse (progression-free survival/PFS) (Fig. 5a, b). The discriminations of OS and PFS are with time-dependent area under the curve (tAUC) value of 94.6%-91.3% in the time frame of 12.4 month (M) to 57.2M (Fig. 5a) and 93.7%-86.7% for 10.8M to 50.6M (Fig. 5b), respectively. Collective evidence supports Overlap66 being a novel and robust prognostic multigene panel for pRCC.
We further validated Overlap66 risk score in stratification of pRCC fatality risk using a recently developed R package: contpointr (https://github.com/thie1e/cutpointr). An optimal cutoff point was obtained with Kernel smoothing model coupled with 1000 bootstrapping. This cutoff point classifies pRCC fatality risk at 0.78 sensitivity and 0.84 specificity or the sum of sensitivity and specificity (sum_sens_spec) value of 1.62 (Fig. 6a). Risk stratifications of out-of-bag bootstrap samples (n=1000) occurred most frequently at sum_sens_spec 1.6 (Fig. 6b), which closely approximates sum_sens_spec 1.62 associated with the optimal cutoff point on the full cohort (Fig. 6a). The fatality risk stratifications of the in-bag samples (n=1000, average 63.2% of full samples) and the out-of-bag samples (n=1000) were at the median sum_sens_spec values of 1.62 and 1.60 respectively. Taken together, these bootstrap analyses reveal a good out-of-sample performance of Overlap66 in classification of pRCC fatality risk, supporting Overlap66’s application in real world. This potential is strengthened by the effectiveness of the risk classification with a range of cutoff points (Fig. 6b, c).
We subsequently optimized Overlap66. As OIP5 expression was at 1.9 folds in ACHN OIP5 tumors compared to ACHN EV tumors (Additional file 13: Table S3), we defined a subgroup of DEGs as those with p.adj < 0.05 and fold ≥ |1.9| in ACHN OIP5 tumors compared to ACHN EV tumors. These DEGs (n=298) share 21 overlap genes (Overlap21) with primary pRCC-derived DEGs (Additional file 7: Fig. S7c, Additional file 16: Table S6b). As expected, Overlap21 is a sub-group of Overlap66 (Table 1). Overlap21 risk scores predict OS possibility under both UV and MV Cox models with comparable efficiency as Overlap66, evident by Hazard ratio (HR) and 95% confident interval (CI) (Additional file 7: Fig. S7b). Similar prediction efficiencies for PFS between Overlap21 and Overlap66 were also observed (Additional file 7: Fig. S7b). Overlap21 effectively stratifies the risk of mortality and PFS; the discriminations possess high tAUC values (Fig. 5a, b). In comparison, Overlap21 seems marginally less effective compared to Overlap66 in the discriminations of OS and PFS (Fig. 5a, b). Nonetheless, the Overlap21-mediated predictions are clearly effective.
The utility of Overlap21 in assessing pRCC fatality risk is further illustrated by its impressive separation of disease-specific survival (DSS) risk (Fig. 7a, c). DSS is more specific compared to OS in addressing factors contributing to cancer-caused deaths. Overlap66 did not perform well in DSS estimation (data not shown), which might be attributable to the small number of events (disease-specific death n=27) in the context of the large number of variables (n=66 in Overlap66). We thus generated Overlap21plus by using Overlap21 as the basis, and the rest of DEGs within Overlap66 were added if they remain risk factors for decreased OS after adjusting age at diagnosis, sex, and T stages (Additional file 16: Table S6s). However, Overlap21plus was not superior to Overlap21 in the estimation of OS and PFS (data not shown). Nonetheless, the risk score of Overlap21plus predicts DSS risk in a comparable efficiency as Overlap21 (Additional file 7: Fig. S7b); its ability to classify DSS possibility was marginally superior to Overlap21 (Fig. 7a-c).
Instead of using time-dependent ROC (receiver-operating characteristic) in evaluating the performance of Overlap66, Overlap21, and Overlap21plus for their prognostic prediction, we further examined their prediction performance using the intact population (i.e. without the time component) by both ROC-AUC and PR-AUC curves. The precision-recall (PR) curve is used to account for the imbalance nature of dataset; the event rates are 14.6% (41/280) for OS, 18.9% (53/280) for PFS, and 9.6% (27/280) for DSS, which are much less than 50%. PR-curve was suggested to evaluate biomarker’s discriminative performance (44). According to both ROC-AUC and PR-AUC curves, Overlap66 predicts OS and PFS possibilities better than Overlap21 (Fig. 5c, d), while Overlap21plus holds a slight edge over Overlap21 in estimating DSS possibility (Fig. 7d, e).
Alterations in immune cell subsets in high-risk pRCC tumors
Tumor-associated immune cells play critical role in tumor initiation and progression (45, 46), suggesting alterations of immune components in Overlap66-stratified high-risk pRCC tumors compared to those of low-risk. To examine this possibility, we profiled all 22 leukocyte subsets in 280 primary pRCC tumors within the TCGA Pancancer dataset using CIBERSORTx (https://cibersortx.stanford.edu/index.php) (47). Significant alterations in several immune cell subsets between high-risk (n=32) and low-risk tumors (n=248) were detected (Fig. 8). Increases in B naïve cells, T follicular helper cells (Tfh), CD4 T memory (activated) cells, and CD8 T (p = 0.075) cells were detected in high-risk local pRCC tumors (Fig. 8a), indicating persistent immune reactions towards tumors; this scenario is not uncommon, evident by the co-existence of ATM-derived tumor surveillance (antioncogenic actions) with oncogenic actions during cancer initiation and progression (48). However, CD8 T cells expressed an upregulation of programmed cell death protein 1 (PDCD1 or PD1) (Fig. 8b), a major mechanism contributing to CD8 T cell exhaustion in cancer (49). Additionally, T regulatory (Treg) cells suppress T cells activation via downregulation of CD80/86 in antigen-presenting dendritic cells (50) and a significant elevation of Treg cells was observed in high-risk pRCC tumors (Fig. 8a). Alterations in M1 and M2 composition in high-risk pRCCs (Fig. 8a) are consistent with the contributions of tumor-associated macrophages in cancer progression (51). Decreases in macrophages M2 in high risk pRCC tumors is supported by a downregulation of β-2-adrenergic receptor (ADRB2) in these tumors (Fig. 8c); the receptor was associated with M2 macrophages (52). Reductions of activated mast cells in high-risk pRCC tumors (Fig. 8a) suggest a downregulation of immune reactions in facilitating pRCC progression. While B naïve cells, CD8 T cells, M2 macrophages, and activated master cells are similarly clustered in both Overlap66 stratified high- and low-risk pRCCs (Additional file 10: Fig. S10), activated CD4 T memory cells, Tfh, Treg, and M1 macrophages in the high-risk tumors display different clustering patterns from their counterparts in the low-risk pRCCs (Fig. 8d-g). Collectively, changes in immune components in high-risk pRCC tumors stratified by Overlap66 risk scores favor the development of an immune suppressive microenvironment, which might be a mechanism underpinning pRCC progression. This concept provides an additional evidence supporting Overlap66 being a novel and effective prognostic biomarker for pRCC.
Critical contributions of PLK1 to OIP5-promoted growth of pRCC tumors
PLK1 (Polo-like kinase 1) is one of the upregulated DEGs identified in relation to OIP5 upregulation in both xenograft tumors and primary pRCC, i.e. a component gene of Overlap66 (Table 1). Its upregulation in xenografts produced by ACHN OIP5 cells compared to those derived from ACHN EV cells was demonstrated by RNA-seq and real-time PCR (Fig. 9a, b). In primary pRCC tumors, OIP5 expression correlates with PLK1 expression with a Pearson correlation value of 0.7 (UALCAN, ualcan.path.uab.edu/home). OIP5, which is also known as Mis18β, is an essential component of the Mis18 complex that is required to load a histone H3 variant CENP-A (centromere protein A) to centromere of newly synthesized DNA strand in early G1 phase (53, 54). PLK1 contributes to CENP-A loading via phosphorylation of M18BP1, a component of the Mis18 complex (55). In line with this knowledge, we examined whether PLK1 kinase activity plays a role in OIP5-promoted pRCC growth.
PLK1 inhibitors have been developed and approved by FDA as Orpha Drug Designation for cancer therapy (56, 57). The PLK1 inhibitor BI2536 caused G2/M arrest with concurrent reduction in G1 phase in ACHN OIP5 cells without apparent effects on cell cycle distributions of ACHN EV cells at the conditions used (Fig. 9c). We then treated mice bearing ACHN EV or ACHN OIP5 cell-produced xenograft tumors when tumors with BI2536 reached 100mm3. In the vehicle treatment group, the OIP5 tumors grew significantly faster compared to the EV tumor (Fig. 9d). Administration of BI2536 had no effects on the growth of ACHN EV tumors but significantly inhibited the growth of ACHN OIP5 tumors (Fig. 9e, f). In the presence of BI2536, ACHN OIP5 tumor showed marginally slower growth compared to ACHN EV tumors (Fig. 9g). Inhibition of PLK1 significantly increases the survival of mice bearing ACHN OIP5 tumor (Fig. 9h). As ACHN is a metastatic pRCC cell line (39), evidence supports inhibition of PLK1 being an option in treating metastatic pRCCs with OIP5 upregulation. Collectively, the above observations indicate synthetic lethality between OIP5 and PLK1 in metastatic pRCCs.