|Home | About | Journals | Submit | Contact Us | Français|
Papillary renal cell carcinoma (PRCC) is the second most common renal cell carcinoma (RCC) that can be further subdivided into type 1 (PRCC1) and type 2 (PRCC2) RCCs based on histological and genetic features. PRCC2 is often more aggressive than PRCC1. While integrin-associated protein complexes mediate tumorigenesis and metastases in many types of cancers it is not known whether integrin-mediated signaling impacts PRCC and differs between PRCC1 and PRCC2. In this study, we combined the analysis of five PRCC gene expression datasets derived from Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) by using integrative bioinformatics pipelines. We found 1475 differentially expressed genes among which 37 genes were associated with integrin pathways. In comparison with PRCC1, PRCC2 cases showed upregulated expression of α5-integrin (ITGA5) whereas the expression of α6- (ITGA6) and β8-integrins (ITGB8) was downregulated. Because PRCC2 occurs more frequently in men, the meta-analysis was extended to explore the gender effects. This analysis revealed 8 genes but none of them was related to integrin pathways suggesting that other mechanisms than integrin-mediated signaling underlie the observed gender differences in the pathogenicity of PRCC2.
The second most common histological subtype of renal cell carcinomas (RCC) is the papillary RCC that accounts for 10–20 percent of all renal cancer cases . Papillary renal cell carcinoma (PRCC) has two subtypes defined by different histological features. PRCC1 shows both papillae and tubular structures covered by small cells with scanty cytoplasm and small oval nuclei. PRCC2 indicates only papillary structures covered by large cells with abundant eosinophilic cytoplasm and large, spherical nuclei with prominent nucleoli . PRCC2 is often a more aggressive disease that is associated with less differentiated histology phenotype, high number of nodal and distant metastases and worse survival rates comparing with PRCC1 [3–4]. PRCC1 and PRCC2 have distinct genetic backgrounds . While overexpression or activating mutations of MET proto-oncogene encoding for a hepatocyte growth factor receptor (HGFR) are common in PRCC1, PRCC2 has been associated with activation of the NRF2-ARE pathway and CDKN2A silencing [6–8]. Interestingly, differential expression of components regulating cell-extracellular matrix (ECM) interactions has also been implicated in PRCC1 . However, how cell-ECM interactions are modulated in PRCC2 remains unknown.
Integrins are heterodimeric cell surface receptors that mediate interactions between cells and the ECM . In solid tumors integrins regulate cancer initiation, stemness, drug resistance and metastasis by regulating the assembly of large multiprotein complexes . These integrin-mediated complexes not only facilitate cells to adhere to the ECM but also convey various signals between cells and their immediate microenvironment to regulate cell behavior [10–11]. Jones et al. reported that abnormal chemokine receptor signaling modulates the activity of α3-, α5-, β1-, β3- and β4-integrins thereby regulating the adhesive properties of clear cell renal cell carcinoma (ccRCC) cells in vitro . However, possible contributions of modified integrin pathways underlying the histological differences between PRCC1 and PRCC2 have not been thoroughly addressed. Intriguingly, PRCC2 occurs more frequently and is more aggressive in male patients when compared with female patients, but the mechanism remains largely unknown [13–15]. In order to address these issues, we combined available RNA-seq and microarray datasets of PRCC samples from GEO and TCGA and compared the mRNA expression profiles of PRCC1 and PRCC2. To get a comprehensive insight into molecular mechanisms and differences between PRCC1 and PRCC2, we performed a meta-analysis of gene expression focusing on the integrin-associated pathways using these large GEO and TCGA datasets. Our analysis highlighted differential regulation of integrin pathways between PRCC1 and PRCC2 but found no significant gender-associated changes in integrin pathway genes between male and female patients with PRCC2.
Four microarray datasets of PRCC patient samples for which matched clinical information was available were obtained from GEO by using GEOquery . In addition, one RNA-seq dataset of PRCC samples was obtained from TCGA-kidney renal papillary cell carcinoma (KIRP) database by using TCGA-assembler  (Figure (Figure1A1A and Table Table1).1). After removing non-PRCC1/PRCC2 (or unidentified) samples and PRCC cases with mixed subtyping, a total of 138 PRCC1 and 135 PRCC2 samples were selected for further analysis (Figure 1A–1C). PCA biplots of the six quality control criteria and the five PRCC datasets supported the inclusion of all of the datasets for meta-analysis (Figure (Figure2A2A and Table Table2).2). Three main meta-analysis methods by combining p-value in MetaDE package were employed: 1) Maximum p-value (maxP), 2) Minimum p-value (minP) and 3) r-th ordered p-value (roP) (Figure (Figure2B)2B) . 1758, 1558 and 1976 differentially expressed (DE) genes were detected by maxP, minP and roP evaluation criteria, respectively, using detection competency curves and false discovery rate (FDR) cut-off less than 0.05 (Figure 2B, 2C and Table Table3).3). All of the three analyses highlighted integrin pathways with significant overlap such that 37% of the DE genes within integrin pathways were shared by all three analyses and 55% were shared between two out of three analyses (Supplementary Figure S1A). Meta-analysis performed with maxP criteria, however, identified more integrin pathway-related DE genes than the other two and was thus selected for further analysis. All the 1758 DE genes obtained with maxP criteria were analyzed by using meta-analysis for pathway enrichment (MAPE) within the MetaPath package to reveal cellular pathways differentially regulated between PRCC1 and PRCC2 . MAPE analysis identified 115 enriched pathways when the MAPE_I method (threshold set to 0.2) which integrates gene and pathway information from the DE genes was employed (Figure (Figure3)3) . Cell adhesion was one of the pathways appearing from this analysis but not among the most prominently enriched pathways. To further refine the list of DE genes between the PRCC subtypes we employed a random-effect model (REM) meta-analysis approach that estimates the differences in gene expression across all the different datasets by combining the individual effects sizes into a meta-effect size (ES) [18–20]. To this end MetaDE.ES was applied for the raw data in the five datasets (Table (Table1)1) resulting in 2610 DE genes when the FDR threshold was set to 0.05 (Supplementary Figure S1B and Table Table3).3). Finally, the DE genes obtained from each of these two meta-analyses were plotted into a Venn diagram that revealed 1475 DE genes common to both maxP meta-analysis method that combines p-values and REM meta-analysis method that combines effect sizes (Figure (Figure4A4A).
The 1475 DE genes obtained from the overlay of p-value and ES combined meta-analysis were further examined using PANTHER gene ontology-slim biological process analysis. This analysis revealed that metabolic process and cellular process were among the top enriched categories (Figure (Figure4B).4B). PANTHER pathway ontology of the 1475 DE genes demonstrated 37 genes related to integrin pathway which was the second most enriched pathway (Figure (Figure4C).4C). These results were corroborated by similarly overlapping the DE genes from each of the two remaining meta-analyses based on p-values (minP and roP) with the DE genes obtained from REM meta-analysis. As expected, PANTHER pathway ontology analysis found integrin signaling pathway to be enriched within shared DE genes between REM and minP (1265 genes) and REM and roP (1300 genes) meta-analyses (Supplementary Figure S1C and S1E).
A functional annotation clustering 2D view report of the integrin pathway associated genes from the maxP/REM overlapping dataset was visualized by using DAVID functional annotation tool (Figure (Figure4D).4D). These genes included α5 (ITGA5)- and β8 (ITGB8)-integrins, both of which are RGD-binding receptors associated with epithelial-to-mesenchymal transition (EMT) [21–22]. Moreover, α6-integrin (ITGA6), a laminin receptor that regulates epithelial cell polarity and growth was found to be differentially regulated between PRCC1 and PRCC2 (Table (Table4)4) [24–25]. Finally, STRING interaction network analysis was performed with these 37 integrin pathway-related genes. This analysis revealed highly connected functional protein-protein interaction network that included all but one of the 37 DE genes (Figure (Figure4E).4E). To relate the differential expression of the three integrin genes in PRCC patients to healthy controls we extracted the expression data from two studies, the TCGA-KIRP study based on RNA-seq data which was also included in our meta-analysis and an independent study by Jones et al. (GSE15641) [7, 26]. This analysis revealed that ITGA5 was significantly downregulated in PRCC1 patients when compared with healthy controls (Figure (Figure5A).5A). In contrast, a tendency for modestly elevated levels of ITGA5 expression was noted for PRCC2 although this was not statistically significant (Figure (Figure5A).5A). ITGA5 expression levels in unclassified PRCC patients in the GSE15641 dataset tended to be higher than in healthy controls but this difference was not statistically significant (Figure (Figure5B).5B). Both studies were consistent with significant downregulation of ITGA6 expression that was particularly evident in PRCC2 patients (Figure (Figure5C).5C). The GSE15641 dataset also displayed a robust downregulation of ITGA6 in the PRCC samples (Figure (Figure5D).5D). In contrast to reduced ITGA6 levels, upregulation of ITGB8 expression was observed in PRCC patients with the highest expression levels seen in PRCC1 patients (Figure (Figure5E).5E). The GSE15641 dataset similarly showed significant upregulation of ITGB8 in PRCC patient samples (Figure (Figure5F).5F). Taken together, the combined meta-analysis extracted three candidate integrin genes whose differential regulation may underlie the different pathogenic properties of the two PRCC subtypes. Comparing with TCGA-KIRP study, the expression pattern of these three integrins predicts that the majority of the cancer samples in the GSE15641 dataset, where matched clinical information was not available, belong to the PRCC subtype 2.
Clinical data shows that PRCC2 is more frequent and tends to be more aggressive in male patients . To address the possible underlying genetic factors we subjected male and female PRCC2 patient data for meta-analysis to study the DE genes. Given that information on gender was not available for two GEO datasets, the remaining three datasets (GSE11024, GSE2748 and TCGA-KIRP) were included for further analysis (Figure (Figure1C).1C). After removing PRCC1 cases the selected datasets consisted of 34 female and 73 male PRCC2 samples. To provide objective quality control and inclusion/exclusion criteria of the filtered datasets for meta-analysis, we performed MetaQC and runQC packages in R. All of the three datasets passed the quality control test and were included for further analysis (Figure (Figure6A6A and Table Table5).5). maxP meta-analysis detected 8 genes that were differentially regulated between female and male PRCC2 samples using a FDR cut-off under 0.05 (Figure (Figure6B).6B). The expression patterns of the 8 genes in the patient samples were visualized in a heat map (Figure (Figure6C).6C). The PANTHER GO-Slim biological process analysis revealed that metabolic process was the major functional category associated with the 8 genes (Figure (Figure6D).6D). STRING Interaction Network analysis of the extracted genes showed that RPS4Y1 interacted with translation elongation factor EEF1A2 (Figure (Figure6E).6E). However, RPS4Y1 is a Y-chromosome-specific gene and when a similar meta-analysis addressing the gender-dependent differences in gene expression was performed for PRCC1 patients where no gender bias has been reported, RPS4Y1 was also detected among the DE genes (Supplementary Figure 2). Therefore, it is unlikely that elevated RPS4Y1 expression is associated with the more aggressive pathogenesis of PRCC2 in male patients. In addition to RSP4Y1, the analysis of gender-specific DE genes in PRCC2 patients also highlighted interaction of two members of the membrane-spanning 4-domain family, subfamily A (MS4A) transmembrane (TMEM)-176A and 176B that have been found to be deregulated in multiple cancer types . However, none of the genes were associated with cell adhesion or integrin pathways.
High-throughput genomic datasets of PRCC patient material have accumulated in the public databases such as GEO and TCGA. Meta-analysis combining these datasets helps to increase the statistical power in mining and explaining the underlying mechanisms driving pathogenesis of PRCC subtypes. Here a meta-analysis was performed to extract the differentially expressed genes between PRCC1 and PRCC2 with a particular focus on the integrin-related pathways. The combined meta-analysis extracted 1475 DE genes, 37 of which were associated with integrin pathway. The more aggressive PRCC2 was associated with modest upregulation of ITGA5, a fibronectin receptor known to promote tumor progression in ccRCC . Cell biological data supports a positive role for α5β1-integrin in promoting mesenchymal cell migration and cancer cell invasion [30–31]. Ligation of α5β1-integrins have also been reported to support cell survival in suboptimal growth conditions [32–33]. It is likely that upregulated ITGA5 expression in part contributes to the metastatic properties of PRCC2. ITGA6 expression was significantly downregulated in PRCC samples, especially in PRCC2, when compared with normal kidneys (Figure (Figure5C5C and and5D).5D). α6-integrin is a laminin receptor which delineates the basal membrane of polarized epithelial cells and synergizes with growth signals to support cell proliferation and migration . Twist-mediated transformation of renal cancer cells led to reduced expression of ITGA6 . Reduced expression of α6β4-integrin may enhance cancer cell dissemination from primary tumors and has been reported to positively correlate with prostate cancer progression . However, increased expression of ITGA6 was reported for metastatic cells of many other solid cancer types . Interestingly, α6-integrin is highly expressed in stem cells including cancer stem cells . Our finding on decreased ITGA6 levels raise the possibility that enhanced dissemination of kidney cancer cells could in part explain the aggressive pathogenesis of PRCC2. However, further functional studies are required to confirm these findings. ITGB8 was the most robustly upregulated integrin in both PRCC datasets. Especially high levels of ITGB8 were seen in PRCC1 samples. ITGB8 forms an αVβ8-integrin heterodimer that plays a critical role in activation of latent TGF-β . TGFβ-activation in turn drives epithelial-to-mesenchymal transition which may directly contribute to cancer cell migration and growth. TGFβ can function as both tumor suppressor and oncogene depending on other coinciding signals . How elevated TGFβ-signaling might regulate PRCC pathogenesis is not known but it has been reported that TGFβ activation modulates metastatic properties of some RCCs [42–43]. TGFβ activation appears to suppress growth of ccRCC tumors but the effect on PRCC remains poorly understood .
A meta-analysis was performed to identify gender-specific genes which may explain the different incidence and aggressiveness of PRCC2 between male and female patients. RPS4Y1, EEF1A2, TMEM176A and TMEM176B were preferentially up-regulated in male patients. RPS4Y1, a Y chromosome-linked gene , was found to functionally interact with a translation factor EEF1A2. While these proteins were linked to metabolic processes EEF1A2 may indirectly associate with cell adhesion related pathways as it has been shown to regulate MMP-9 expression and thereby influence migration and metastatic properties of pancreatic cancer cells . Our meta-analysis approach revealed three key integrin subunits that were differentially associated with PRCC subtypes. An important limitation of the current study is that integrin activation and stability is frequently regulated at post-translational level. Therefore inclusion of proteomic data as well as functional cell biological studies are required for future studies to confirm and validate the specific roles of α5-, α6- and β8-integrins highlighted here by genetic meta-analysis.
173 datasets were found in GEO series (http://www.ncbi.nlm.nih.gov/gds) by search term “papillary renal cell carcinoma” (Searched on 10th of February 2016). After manually removing datasets in which PRCC1 and PRCC2 had not been separately classified, four datasets were obtained (Table (Table1).1). For the GSE26574 dataset, 22 PRCC1 and 12 PRCC2 samples were analyzed by expression profiling microarray using a GPL11433 platform without gender information . In GSE7023 dataset, 14 PRCC1 and 16 PRCC2 samples were analyzed by expression profiling microarray using GPL4866 platform without gender information . In GSE11024 and GSE2748 datasets, gender information was available for a total of 25 PRCC1 and 21 PRCC2 samples which were analyzed by expression profiling microarray using GPL6671 and GPL570 platforms, respectively [49–50]. In addition to microarrays datasets, one large PRCC dataset was available via TCGA-KIRP and consisted 77 PRCC1 and 86 PRCC2 with gender and clinical information . This dataset had been analyzed by RNA-seq using Illumina Hiseq RNASeqV2-platform (http://cancergenome.nih.gov/).
Capacity computing environment (Finland CSC Taito-shell application server: https://research.csc.fi) was used for running integrative bioinformatics pipelines. For the four GEO datasets, the raw gene expression data and clinical information was downloaded from GEO by using GEOquery package in R (BioConductor; https://www.bioconductor.org/). GEOquery is an open-source and open-development software project [52, 16]. TCGA-assembler was applied for downloading and processing the TCGA-KIRP gene expression raw data and clinical information [53, 17]. After removing non-PRCC samples from the datasets, the remaining data was saved to a comma-delimited text-file by using Excel as described previously . Objective quality control analysis was performed by importing the processed raw datasets in log-format into R by using MetaDE . Largest interquartile range (IQR) method was used to calculate gene-wise expression. Non-expressed (30%) and non-informative (30%) genes were filtered out. MetaQC was applied to include or exclude the processed datasets as described previously [55, 18, 56].
The MetaQC processed datasets were subjected to meta-analysis by combining p-values and effect sizes . For the meta-analysis, the five datasets were combined and subjected to quality control and inclusion/exclusion criteria. These criteria consisted of internal quality control (IQC) index evaluating the homogeneity of coexpression and external quality control (EQC) index supervised by external pathway information, accuracy quality control indexes for genes (AQCg) and pathways (AQCp) and consistency of differential expression quality control (CQCg and CQCp) indexes that collectively depict the reproducibility and consistency of the data between the different individual the studies and results from the combined meta-analysis . Three different meta-analysis methods for combining p-value; maxP, minP and roP were applied resulting in 1558, 1976 and 1758 differentially expressed (DE) genes, respectively (FDR = 0.05, Table Table3)3) [20, 18]. For visualization, a heatmap plot of the DE genes (FDR = 0.05) was created by using MetaDE. A heatmap for pathway enrichment was provided by MetaPath packages, and MAPE_I was under q-value = 0.2 threshold . Additionally, 2610 genes were obtained from a “REM” meta-analysis of the same raw data by combining effect sizes using MetaDE (FDR = 0.05) . A Venn diagram was created to display an overlay of 1475 DE genes using Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/) . The 1475 DE genes extracted from the overlay meta-analysis were further analyzed with DAVID 6.7 and PANTHER gene ontology to classify gene pathway, functions and interactions (https://david.ncifcrf.gov/) [60–61]. Protein-protein interactions of the identified DE genes were also analyzed using STRING v10 online tool that visualizes known and predicted protein-protein interactions (http://string.embl.de/) [62–63]. The gender comparison of PRCC2 samples were processed with MetaQC and MetaDE by combining the p-values of the raw data after PRCC1 cases were manually removed. PANTHER GO.Slim Molecular function analyses were performed as described earlier (http://pantherdb.org/) . Expression analysis of ITGA5, ITGA6 and ITGB8 in Renal Cell Carcinomas and normal Kidney was obtained from Oncomine  and TCGA-KIRP.
We thank Dr. Xingbin Wang and Dr. Dongwan D. Kang, Department of Human Genetics, University of Pittsburgh, for suggestions and instructions with the MetaOmics R packages.
CONFLICTS OF INTEREST
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.
This work was funded by Academy of Finland (140974, 263770, and 135560) and Biocenter Oulu.