Genome-wide pro ling identi cation of prognostic novel signatures in papillary renal cell carcinoma based on large-scale sequencing data

Xi Tian Fudan University Shanghai Cancer Center https://orcid.org/0000-0003-1965-0647 Wen-Hao Xu Fudan University Shanghai Cancer Center Aihetaimujiang anwaier Fudan University Shanghai Cancer Center Hong-Kai Wang Fudan University Shanghai Cancer Center Fang-Ning Wan Fudan University Shanghai Cancer Center Yu Zhu Fudan University Shanghai Cancer Center Da-Long Cao Fudan University Shanghai Cancer Center Yi-Ping Zhu Fudan University Shanghai Cancer Center Guo-Hai Shi Fudan University Shanghai Cancer Center Yuan-Yuan Qu (  quyy1987@163.com ) https://orcid.org/0000-0003-1179-2800 Hai-Liang Zhang Fudan University Shanghai Cancer Center Ding-Wei Ye Fudan University Shanghai Cancer Center


Background
Papillary renal cell carcinoma (PRCC), the second most common histologic subtype of renal cell carcinoma, originates from tubular epithelial cells, accounts for about 10%-15% of all renal tumors [1,2].
The histological features of PRCC are brovascular cores with papillary arrangement of tumor cells [2]. PRCC can be divided into two subtypes: type I and type II. Type 1 PRCC was characterized by monolayer, light staining and basophilic small cells, often with rich foam macrophage in ltration, and type II showed high Fuhrman grade, eosinophilic cytoplasm and pseudostrati ed arrangement in the center of the nipple [3].
Recent studies have shown that type II PRCC is a heterogeneous tumor, which can be further subdivided into other subtypes according to the genetic and molecular composition of the tumor, re ecting different clinical course and prognosis [4]. It is reported that PRCC type II patients have worse outcome than clear cell renal cell carcinoma (ccRCC) [5,6] and it is di cult to predict the prognosis of patients. Some genetic features of PRCC have been recognized such as mutations of MET, NF2, SETD2, and Nrf2 pathway genes [7][8][9]. However, these mutations were found in only 10-15% of PRCC tumors in these studies [10,11]. In recent years, with better understanding of RCC molecular biology, targeted drugs and therapeutic effects have been improved [12].
With the rapid development of microarray technology and bioinformatics analysis, our understanding of differentially expressed genes (DEGs) and functional pathways related to the occurrence and development of PRCC has become more comprehensive. However, the rarity of this tumor has become an obstacle in identifying potential markers to distinguish PRCC and provide potential therapeutic targets. In this study, two mRNA microarray data sets were downloaded from Gene Expression Omnibus (GEO) [13] for analysis to obtain DEGs between cancer tissues and normal adjacent tissues. Subsequently, in order to further understand the molecular mechanism of tumorigenesis, gene ontology (GO) [14] and Kyoto genome encyclopedia (KEGG) [15] pathway enrichment analysis were carried out. The protein-protein interaction (PPI) network reveals the speci c functions of all proteins and describes the importance of these interactions in biological processes, molecular functions, and signal transduction.
To determine the candidate biomarkers and their possible role in PRCC, this work focused on analyzing gene expression pro les, reveal potential biointeraction networks, and assessing prognostic value. We speculated that the oncogenic activity of important hub genes was associated with poor prognosis, which may be a potential therapeutic target for PRCC.

Materials And Methods
Original biological microarray data Gene Expression Omnibus (GEO) is a public functional genomic database that stores high throughout gene expression data, chips and microarrays. The Original DNA microarray data were obtained from GEO, for patients with PRCC. Two chip datasets GSE48352, GSE26574 were downloaded from GEO (Affymetrix GPL16311 platform, Affymetrix GPL11433 platform, respectively). The corresponding genes transformed into a probe were converted into a symbol according to the annotation information on the platform.

Screening and identi cation of DEGs
The DEGs between PRCC and non-cancerous samples were screened and identi ed by GEO2R.
Delineating parameters such as P-values, Benjamini and Hochberg false discovery rate (FDR) and fold change were used to lter DEGs and applied to provide a balance between the discovery of statistically important genes and the limitations of false positives. A probe set that does not have a corresponding gene symbol or a gene that has multiple probe sets is knocked out or averaged. logFC (fold change) > = 1 or <=-1 and P-value < 0.01 were considered statistically signi cant.

PPI network construction and module analysis
In this study, we used Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 10.0) online database to predict PPI networks of DEGs and analyze functional interactions between proteins [16]. This may help to further understand the underlying mechanisms of the development and progression of PRCC. Cytoscape (version 3.5), an open source bioinformatics software platform, was designed to visualize molecular interaction networks [17]. Cytoscape's plug-in Molecular Complex Detection (MCODE) (version 1.4.2) can cluster a given network based on topology to nd dense connection areas [18]. The most important modules in the PPI network are selected as follows: MCODE Score > 24. Subsequently, KEGG and GO analysis of the genes in the module was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov; Version 6.8) [19].

Hub genes selection and analysis
After applying the MCODE plug-in, a network of the 38 genes and their co-expression genes was analyzed using cBioPortal (http://www.cbioportal.org) online platform [20]. ClueGO, a Cytoscape plug-in that could visualize the non-redundant biological terms for large clusters of genes in a functionally grouped network [21]. The biological processes from GO and KEGG pathway analysis of hub genes were visualized  [22,23]. Heat map, based on hierarchical clustering algorithm, was constructed using the phenotype and gene expression pro le of 323 samples in TCGA.

Functional enrichment of DEGs
Biological properties such as biological processes (BP), molecular functions (MF), and cellular components (CC) were extracted from gene ontology (GO) enrichment analysis to determine the role of DEGs in PRCC. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies. DAVID was applied to discern the role of developmentrelated signaling pathways in PRCC. P-values < 0.05 are considered statistically signi cant. GO and KEGG enrichment was analyzed and displayed using bubble chart.

Statistical analysis of TCGA cohort
Phenotype and expression pro les of hub genes in 323 PRCC patients from TCGA were analyzed and displayed to predict prognostic value. The expression levels of each gene were divided into three groups: high, medium and low. Subsequently, to compare the outcome of high group and low group, the univariate survival analysis of the 8 hub genes was performed using Kaplan-Meier curve. Univariate and multivariate analyses were performed with Cox logistic regression models to nd independent variables, including age (ref. values were taken using X-tile software [24].

Immunohistochemical (IHC) staining and evaluation
Rabbit anti-PTTG1 monoclonal antibody was used (ab128040, Abcam, USA). The positive or negative staining of a certain protein in a FFPE slide was independently evaluated by two experienced pathologists and determined as follows. The overall IHC score from 0 to 12 was evaluated according to the multiply of the staining intensity and extent score, as previously described [25].
Statistical analysis of FUSCC cohort and potential networks of the hub genes According to the IHC score, the patients were divided into two groups: high expression group and low expression group of PTTG1. Correlation analyses between the expression of PTTG1 and clinicopathological features were carried out by chi-square test. In order to explore the factors related to

Results
This study is divided into three phases. In the rst phase, we used the information in GEO platform to evaluate DEGs; in the second phase, we constructed a PPI network and evaluated their interactions speci city based on co-expression and functional annotation. In the third phase, phenotype and expression pro les of hub genes in 323 PRCC patients from TCGA were analyzed and displayed to predict prognostic value.

Identi cation of DEGs in PRCC
After normalizing the microarray results, we identi ed DEGs, including 1,270 probes in GSE48352, 826 probes in GSE26574. As shown in the Venn diagram, the overlap in both two datasets contained 473 differentially expressed genes between tumor tissues and adjacent normal tissues (Fig. 1A).
PPI network construction and module analysis PPI network of DEGs was built (Fig. 1B) and the most important modules were identi ed using Cytoscape plugin (Fig. 1C). Functional analysis of the 38 genes involved in the module was performed using DAVID. The functional analyses of 38 genes involved in this module were analyzed using DAVID. Enrichment pro les showed that hub genes in this module were primarily enriched in sister chromatid segregation (49.37%), mitotic nuclear division (13.92%), anaphase-promoting complex-dependent catabolic process (11.39%; Supplementary Fig. 1). Through functional analysis, we subjectively selected eight genes related to the cell cycle and highlighted them in the PPI network of 38 hub genes.  Table 1, modules with both P value < 0.05 and FDR < 0.1 were considered signi cant.  (Fig. 4A). Univariate Cox regression analyses was carried out, expression levels of the eight hub genes and traditional prognostic factors such as AJCC stage and pTNM stage were considered signi cantly related to OS (p < 0.05; Fig. 4B). In multivariate Cox regression analysis models, traditional prognostic factors were not always relevant to OS. Expression levels of the eight hub genes were still signi cantly correlated with OS (p < 0.05; Fig. 5A, C, E, G, I, K, M, O) and the AUC indices for the TCGA-OS were from 0.686 to 0.819 (Fig. 5B, D, F, H, J, L, N, P). After integrating all the signi cant clinicopathological parameters and gene expression pro les in the Cox regression models, we generated eight formulas for OS (Supplementary table 1).

PTTG1 expression and clinicopathological characteristics in XXX PRCC patients from the FUSCC cohort
Although the AUC of MCM5 is the highest in the eight genes, we did not choose to verify the prognostic value of MCM5 due to the little correlation with PFS. As the AUC of PTTG1 ranked in second place and it correlated signi cantly with both OS and PFS, we chose it to make further veri cation. In order to evaluate the expression level of PTTG1 protein in FUSCC tumor samples, we performed IHC staining and found that the density and intensity of PTTG1 expression in PRCC tissues in the FUSCC cohort were signi cantly higher than those in adjacent normal renal tissues (Fig. 6A-C).  (Table 2). *P value less than 0.05 was considered as statistically signi cant and marked in bold.

Cox regression analyses of FUSCC cohorts and potential networks of the hub genes
In univariate Cox regression analysis models, pTNM stage, AJCC stage, and Furhman grade, tumor size, BMI, were signi cantly relevant to PFS (p < 0.05; Table 3) and OS (p < 0.05; Table 3). Importantly, PTTG1 expression markedly correlated with poor PFS (HR = 2.46, p < 0.01) and poor OS (HR = 2.90, p < 0.01). The potential networks associated with functional annotations of the hub genes include inhibition, protein interaction, ubiquitination, phosphorylation, activating expression, activation, indirect relation, dephosphorylation (Fig. 7). *P value less than 0.05 was considered as statistically signi cant and marked in bold.

Discussion
The concept of PRCC was rst proposed by Mancilla-Jimenez in 1976, thirty-four cases of RCC showed papillary structures. Of these, 85.3% PRCC patient have a better prognosis than other types of RCC [26].
Since PRCC is relatively rare in clinical practice and has been rarely studied, major molecular mechanisms in the pathogenesis are poorly understood. Therefore, potential biomarkers for e cient diagnosis and treatment are urgently needed. In this study A total of 473 DEGs and 38 hub genes were identi ed by microarray data analysis. Among the 38 hub genes, 8 genes relating to cell cycle including BUB1B, CCNB1, CCNB2, MAD2L1, TTK, CDC20, PTTG1, MCM5 were subjectively selected. After statistical analysis, the 8 genes showed clear prognostic value.
BUB1B (spindle detection point protein, also known as BUBR1) is an important functional protein at the detection point of mitosis and the change of BUB1B expression plays an important role in tumorigenesis and progression [27]. Studies have found that BUB1B is overexpressed in kidney cancer and breast cancer, its mutation and overexpression is strongly correlated with Chromosomal instability [27,28]. Yet the prognostic value of BUB1B in PRCC has been rarely reported.
As a member of the cell cycle family, CCNB1 is one of the key factors related to cell detection points [29,30]. Currently, Cyclin B1 overexpression has been found in a variety of human tumors, such as esophageal cancer, non-small cell lung cancer, tongue cancer, and is related to tumor grade, differentiation, invasion and metastasis and prognosis [31]. Thus, there is enough evidence to doubt the role of CCNB1 in human PRCC as an oncogene.
It has been reported that CCNB2 is highly expressed in tumor tissues, such as breast cancer [32], adrenal cortical carcinoma [33], colorectal adenocarcinoma [34]and pituitary adenoma [35]. It has also been reported that serum circulating CCNB2 mRNA level in cancer patients is signi cantly higher than that in normal population and benign diseases [36]. Thus, CCNB2 may have play an important role in the generation and development of PRCC.
It has been shown that interrupting the function of MAD2L1 in mammalian cells can affect the process of spindle examination and lead to the development of aneuploid cells or tumors. The deletion of MAD2L1 gene can cause chromosome instability and drive the development of tumor. In mouse models, deletion of MAD2L1 gene can cause liver and lung cancer [37].
TTK is the basic component of spindle assembly checkpoint(SAC), it plays an important role in the replication of mitotic centrosomes and the correct separation of chromosomes [38]. In order to maintain the division and proliferation of tumor cells, TTK was highly expressed in tumor cells to maintain the normal function of SAC. After inhibiting the function of TTK, SAC is damaged, errors in mitotic metaphase cannot be detected, chromosomes cannot be separated into daughter cells on average, and heteroploidy is further increased, exceeding a certain threshold will cause tumor cell apoptosis, so TTK can serve as an effective anti-tumor target [39,40].
Multiple studies have shown that CDC20 could degrade several important substrate factors to regulate cell cycle progression including Securin [41], Cyclin A [42,43], p21 [44] and Mcl-1 [45]. The protein p21 is considered to be an effector of various tumor inhibition signaling pathways, partly deactivate Cyclin dependent kinases to promote anti -tumor proliferation. Downregulation of p21 expression was detected in a variety of human malignancies, so CDC20 may play its carcinogenic role in part by degrading tumor suppressor protein p21. Thus, the generation of PRCC may be caused part by CDC20.
MCMs protein is closely related to the cell cycle. It is a promoter of DNA replication and plays a key role in regulating cells going from G0 phase to S phase [46]. As one of the MCM protein families, MCM5 protein is closely related to cell proliferation. The expression of MCMs protein can be indirectly understood through the detection of MCM5. Williams et al. reported the application of MCM5 protein antibody to diagnose abnormal prodromal malignant cells in pap cervical smear, and found that the MCM had high sensitivity and speci city in detecting predromal malignant cells in cervical smear [47]. Going has proposed this conclusion in esophageal tissues [48]. These results indicate that MCM5 protein has some predictive value.
PTTG1 is a tumor transforming gene, which can cause cell transformation without the participation of any auxiliary gene, and is closely related to the occurrence of many tumors [49]. PTTG1 has been identi ed as an oncogene. The expression level of PTTG1 is closely related to tumor formation, angiogenesis and metastasis [50]. In our study, we detected the expression level of PTTG1 protein in PRCC from FUSCC cohort by using immunohistochemical staining and we found that PTTG1 protein increased in tumor tissues, and higher expression of PTTG1 is signi cantly relevant to both OS and PFS.
There are still several limitations in this study. First, the data utilized in the study consisted of unbalanced PRCC and control normal samples, which were restricted in quantity and downloaded from the Gene Expression Omnibus database, not generated by us. The chip data contains relatively small PRCC samples in public database and only 323 patients were enrolled form TCGA cohort with corresponding RNA sequence. Second, we only veri ed the expression and prognostic values of PTTG1 in FUSCC cohort and the potential mechanism of signaling pathway in PRCC is not clear, while a series of functional annotations and enrichment analyses have been carried out. The detailed mechanism between the eight genes and PRCC need to be further studied.

Conclusion
In conclusion, DEGs and hub genes may be involved in poor prognosis and recurrence of PRCC in this study. The transcription pro les of BUB1B, CCNB1, CCNB2, MAD2L1, TTK, CDC20, PTTG1 and MCM5 are predictive and may contribute to a better understanding of the potential carcinogenesis and progression of PRCC. Further research needs to elucidate the molecular mechanism and signaling pathway changes of these genes in PRCC.