Abstract
Purpose. Papillary renal cell carcinoma (pRCC) is the second most common histological subtype of adult kidney tumors, with a poor prognosis due to limited understanding of the disease mechanism. Herein, we have performed high-throughput bioinformatic screening to explore and identify potential biomarkers of DNA damage and oxidative stress for pRCC. Methods. RNA sequencing data related to pRCC were downloaded from the TCGA database, and differentially expressed genes (DEG) were identified by a wide variety of clustering and classification algorithms, including self-organized maps (SOM), artificial neural networks (ANN), support vector machines (SVM), fuzzy logic, and hyphenated techniques such as neuro-fuzzy networks. Then DAVID and STRING online biological information tools were used to analyze functional enrichment of the regulatory networks of DEG and construct a protein-protein interaction (PPI) network, and then the Cytoscape software was used to identify hub genes. The importance of key genes was assessed by the analysis of the Kaplan–Meier survival curves using the R software. Lastly, we have analyzed the expression of hub genes of DNA damage and oxidative stress (BDKRB1, NMUR2, PMCH, and SAA1) in pRCC tissues and adjacent normal tissues, as well as the relationship between the expression of hub genes in pRCC tissues and pathological characteristics and prognosis of pRCC patients. Results. A total of 1,992 DEGs for pRCC were identified, with 1,142 upregulated ones and 850 downregulated ones. The DEGs were significantly enriched in activities including DNA damage and oxidative stress, chemical synaptic transmission, an integral component of the membrane, calcium ion binding, and neuroactive ligand-receptor interaction. cytoHubba in the Cytoscape software was used to determine the top 10 hub genes in the PPI network as BDKRB2, NMUR2, NMU, BDKRB1, LPAR5, KNG1, LPAR3, SAA1, MCHR1, PMCH, and NCAPH. Furthermore, the expression level of hub genes BDKRB1, NMUR2, PMCH, and SAA1 in pRCC tissues was significantly higher than that in the adjacent normal tissues. Meanwhile, the expression level of hub genes BDKRB1, NMUR2, PMCH, and SAA1 in pRCC tissues was significantly positively correlated with tumor stage, lymph node metastasis, and the histopathology grade of pRCC. In addition, high expression levels of hub genes BDKRB1, NMUR2, PMCH, and SAA1 were associated with a poor prognosis for patients with pRCC. Univariate and multivariate analyses showed that the expression of hub genes BDKRB1, NMUR2, PMCH, and SAA1 were independent risk factors for the prognosis of patients with pRCC. Conclusion. The results of this analysis suggested that BDKRB1, NMUR2, PMCH, and SAA1 might be potential prognostic biomarkers and novel therapeutic targets for pRCC.
1. Introduction
Renal cell carcinoma (RCC), also known as kidney cancer, is derived from renal tubular epithelial cells and is the most common solid tumor of the kidney, accounting for 3% of adult malignant tumors [1]. It is a heterogeneous group of cancers arising from renal tubular epithelial cells that encompasses 85% of all primary renal neoplasms. Papillary renal cell carcinoma (pRCC) is the second most common histological subtype after clear cell renal cell carcinoma (ccRCC), and 10–15% of RCC histological types are papillary renal cell carcinoma [2]. There are two subtypes of pRCC, type I (basophilic) and type II (acidophilic), and type I has a better prognosis than type II [3]. Most research studies on kidney cancer has focused on ccRCC, and the related studies have shown that compared with ccRCC patients, pRCC patients typically have a lower stage and grade of tumor as well as longer overall survival [4]. The molecular mechanism of pRCC has not been clearly defined. With poor sensitivity to radiotherapy and chemotherapy, surgery is the preferred method for treatment of pRCC, but some patients are prone to metastasis and relapse after surgery. With continued advances in molecular medicine in recent years, the study of the occurrence, development, and metastasis mechanisms of pRCC can help to guide clinical diagnosis and treatment.
The cancer genome atlas (TCGA) project is a joint project of the National Cancer Institute and the National Human Genome Research Institute and aims to apply high-throughput genome analysis technology and to improve the ability to prevent, diagnose, and treat cancer. The cancer genome atlas (TCGA) research network includes analysis of a large number of human tumors to discover molecular aberrations at the DNA, RNA, protein, and epigenetic levels [5]. In this study, TCGA data were used to investigate genes that are deferentially expressed in pRCC. To mine the key genes related to pRCC occurrence and development, we conducted differential gene enrichment (Gene Ontology, GO) analysis and KEGG pathway enrichment analysis, constructed PPI interaction networks, screened hub genes, and performed survival analysis.
2. Materials and Methods
2.1. Data Collection
The published transcriptome data related to papillary renal cell carcinoma were downloaded from TCGA (https://cancergenome.nih.gov/). The data included 289 papillary renal cell carcinoma samples and 32 normal kidney tissues.
2.2. Identification of DEGs
We have performed the edgeR software package in R language (version 3.5.3, https://www.r-project.org/) and a wide variety of clustering and classification algorithms, including self-organized maps (SOM), artificial neural networks (ANN), support vector machines (SVM), fuzzy logic, and hyphenated techniques such as neuro-fuzzy networks to standardize the data and analyze differential expression. Genes with |logFC| > 2.0 and FDR <0.05 were considered differentially expressed genes. To visualize the data graphically, the ggplot2 software package was used.
2.3. GO and KEGG Pathway Analysis
The DAVID database (DAVID; https://david.ncifcrf.gov) was used to perform annotation, visualization, and integrated discovery on the genes identified as significantly differently expressed [6]. Using DAVID, GO analysis was performed, including the analysis of cellular components (CC), molecular functions (MF), and biological process (BP) terms. A value of was considered statistically significant. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.genome.jp/kegg/) is a knowledge base for systematic analysis of gene functions, linking genomic information with higher order functional information [7]. An adjusted value <0.05 was considered statistically significant.
2.4. Hub Genes Selection and Analysis of Modules from PPI Networks
The STRING database (http://string-db.org) aims to provide a critical assessment and integration of protein-protein (PPI) interactions [8]. STRING was used to analyze the selected differentially expressed genes and construct a PPI network. Then, cytoHubba in Cytoscape software (version 3.7.2) was used to screen the top 10 hub genes in the PPI network [9].
2.5. Survival Analyses of Hub Genes
The expression profiles and clinical data of 289 pRCC samples were downloaded from TCGA (http://tcga-data.nci.nih.gov) for the survival analysis of hub genes. The Kaplan–Meier method was used for the survival analysis, and log-rank values were calculated. A log-rank value <0.05 was considered statistically significant.
2.6. Clinical Specimens
A total of 60 paired pRCC samples and adjacent normal renal specimens were collected from Zhuzhou Central Hospital between June 2016 and June 2021. Inclusion criteria for specimen collection: (1) Postoperative pathology examination confirmed pRCC; (2) the patients with neither radiotherapy nor chemotherapy; (3) complete follow-up data were available; (4) the patients understood the purpose and requirements of the study, agreed to participate in the study, and signed a written informed consent, which was reviewed and approved by the Ethics Committee of Zhuzhou Central Hospital.
2.7. Total RNA Isolation and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
The RNA was isolated by TRIzol® reagent (Ambion; USA) from pRCC tissues according to the manufacturer’s protocols. And cDNA was reversely transcribed by PrimeScript RT reagent kit (Takara, China). We conducted RT-qPCR on an ABI 7500 RT-PCR system using the SYBR Premix Ex TaqII Kit (Takara, China). All quantifications were normalized to the level of glyceraldehyde phosphate dehydrogenase (GAPDH) in the reaction. rimers of BDKRB1 was Forward (5′–3′) CAC-TGT-CCT-ACC-GTC-TTT-GTCT, Reverse (5′–3′) CGC-AAA-TCT-TGG-TAG-GTG-GT; NMUR2 forward (5′–3′) GGC-AAG-GCC-ATG-TGT-AAG-ATC, Reverse (5′–3′) GTA-AAA-CGA-CGG-CCAG; PMCH forward (5′–3′) CAC-TGT-CCT-GAC-CGT-CTT-TGT-CT, Reverse (5′–3′) CCA-TAT-GCC-TGT-GGA-GTG-GAA; SAA1 forward (5′–3′) ACC-TGA-GGA-GCC-CCA, Reverse (5′–3′) TCT-GCT-CCT-GGC-AGG-CC.
The comparative threshold cycle (CT) method, which compares the differences in CT values between common reference RNA and target gene RNA, was used to obtain the relative fold changes in gene expression. The expressions were calculated by 2−ΔΔct method. Each experiment was performed in triplicate and repeated three times.
2.8. Statistical Analysis
SPSS 24.0 software was used for statistical analysis, and GraphPad Prism 7.0 software was used for analysis and mapping. All measurement data in the form of mean ± standard deviation (SD), according to two groups and multiple groups of measuring data comparison using Student’s t-tests and one-way ANOVA. The relationship between the RNA expression levels of hub genes BDKRB1, NMUR2, PMCH, and SAA1 in the patients with pRCC tissue samples and the clinical pathological characteristics of patients with pRCC was analyzed through Pearson’s Chi-squared test, and the relationship between the expression of hub genes BDKRB1, NMUR2, PMCH, and SAA1 and the prognosis of pRCC patients was analyzed by Kaplan–Meier survival analysis and the Cox proportional hazard model. was considered to be significantly different.
3. Results
3.1. Identification of DEGs
The data for 289 cases of papillary renal cell carcinoma and 32 cases of normal kidney tissue were downloaded from TCGA and used for this study. The data were normalized and logarithmized, probes without corresponding gene annotation information were removed, and repeated probes were removed to finally get the expression profiles of 17,894 genes and 321 samples. Using the edgeR software package, with |logFC|> 2.0 and FDR <0.05 as the screening conditions for differentially expressed genes, a total of 1,992 DEGs were screened for pRCC, including 1,142 upregulated genes and 850 downregulated genes. Using these selected genes, a volcano map (Figure 1) was generated, and the top 50 gene heat maps with the most significant differences were selected (Figure 1(b)).
(a)
(b)
3.2. GO Term and KEGG Pathway Analyses
In order to better understand the relationships between DEGs and pRCC, we input all DEGs into the online tool DAVID to perform GO analysis. The results revealed that, for GO BP analysis, the DEGs of pRCC were mainly enriched in excretion, epidermis development, ion transmembrane transport, chemical synaptic transmission, chloride transmembrane transport, ion transport, and potassium ion transmembrane transport. For GO CC analysis, DEGs were mainly enriched in integral component of plasma membrane, extracellular region, extracellular space, plasma membrane, apical plasma membrane, anchored component of membrane, proteinaceous extracellular matrix, integral component of membrane, and basolateral plasma membrane. For GO analysis, DEGs were mainly enriched in calcium ion binding, heparin binding, sequence-specific DNA binding, transporter activity, and carbohydrate binding. The GO analysis findings are shown in Figure 2 and Table 1.
(a)
(b)
(c)
We next performed KEGG pathway analysis to analyze the pathways at the functional level. The results showed that DEGs were mainly enriched in neuroactive ligand-receptor interaction, calcium signaling pathway, gastric acid secretion, bile secretion, and pancreatic secretion. The KEGG pathways associated with enriched DEGs associated with pRCC are presented in Figure 2(b) and Table 2.
3.3. Identification of Hub Genes and Analysis of Modules from PPI Networks
The STRING database was used to construct PPI networks for DEGs related to the pathogenesis of papillary renal cell carcinoma. We used the MCODE in Cytoscape software to obtain the main PPI network (Figure 2(c)), and then used cytoHubba in Cytoscape software to identify the top 10 hub genes in the PPI network (Figure 2(c)): recombinant bradykinin receptor B2 (BDKRB2), neuromodulin U receptor 2 (NMUR2), neuromodulin U (NMU), recombinant bradykinin receptor B1 (BDKRB1), lysophosphatidic acid receptor 5 (LPAR5), Kininogen-1 (KNG1), lysophosphatidic acid receptor 3(LPAR3), serum amyloid A1 (SAA1), melanin-concentrating hormone receptor 1 (MCHR1), and precursor melanin-concentrating hormone (PMCH). These 10 hub genes are presented in Figure 2(c).
3.4. Survival Analysis of Hub Genes
Expression data for a total of 289 pRCC samples were downloaded from TCGA. The 10 hub genes were grouped by expression levels, and the data were used to conduct survival analyses. Increased expression levels of BDKRB1, NMUR2, PMCH, and SAA1 were associated with a worse survival rate for pRCC patients (Figure 3).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
3.5. The Expression of Hub Genes BDKRB1, NMUR2, PMCH, and SAA1 in pRCC Tissues and Adjacent Normal Tissues of pRCC Patients
We selected 120 tissue samples (including 60 pRCC tissues and 60 normal adjacent tissues) to analyze the expression of hub genes BDKRB1, NMUR2, PMCH, and SAA1 in pRCC tissues by qRT-PCR. The results showed that the expression of hub genes BDKRB1, NMUR2, PMCH, and SAA1 in pRCC tissues was significantly higher than that in the normal adjacent tissues (Figures 4(a), 4(c), 4(e), and 4(g)). To further investigate the correlation between hub genes BDKRB1, NMUR2, PMCH, and SAA1 expression and pathological features of pRCC, the above samples were divided into high (above the mean) and low (below the mean) hub genes expression groups. Subsequently, the Chi-square test was used to analyze the relationship between hub genes BDKRB1, NMUR2, PMCH, and SAA1 expression level and pathological characteristics of pRCC patients, and the results showed that the expression level of hub genes BDKRB1, NMUR2, PMCH, and SAA1 expression in pRCC tissues were significantly positively correlated with tumor stage, lymph node metastasis, and histopathological grade of pRCC patients (Figures 4(b), 4(d), 4(f), and 4(h)), while the relationship with gender and age of patients was not statistically significant (Tables 3–6).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
3.6. Relationship between Hub Genes BDKRB1, NMUR2, PMCH, and SAA1 Expression and Prognosis of Patients with pRCC
The Kaplan–Meier survival analysis was used to study the relationship between hub genes BDKRB1, NMUR2, PMCH, and SAA1 expression and prognosis of patients with pRCC. The results showed that the overall survival rate of patients with high hub genes BDKRB1, NMUR2, PMCH, and SAA1 expression was significantly lower than that of patients with low hub genes BDKRB1, NMUR2, PMCH, and SAA1 expression (Figure 5). Then we conducted the COX proportional risk model analysis. The univariate and multivariate analyses showed that the expression of hub genes BDKRB1, NMUR2, PMCH, and SAA1 were independent risk factor for prognosis in patients with pRCC (Tables 7–10).
(a)
(b)
(c)
(d)
4. Discussion
Most patients with pRCC have no obvious symptoms or signs at the time of diagnosis, but the disease is often found by B-ultrasound or CT examination during a physical examination. Very few patients exhibit the typical triad signs of kidney cancer: hematuria, abdominal mass, and lumbar pain, and the patients that do exhibit these signs typically have advanced disease. The overall prognosis of pRCC is better than that of ccRCC, but pRCC prognosis is significantly worse than that of ccRCC when pRCC invades the renal vein and/or the inferior vena cava [10]. There is currently no specific treatment for pRCC, and surgical treatment is the first choice in clinical practice. The prognosis of advanced patients is poor, a pRCC is insensitive to radiotherapy and chemotherapy. Therefore, the study of the mechanisms of pRCC development and metastasis will help improve clinical diagnosis and treatment.
In this study, bioinformatics technology was used to mine pRCC transcriptomic data downloaded from TCGA. A total of 1,992 DEGs were identified, including 1,142 upregulated genes and 850 downregulated genes. We performed GO and KEGG pathway enrichment analyses to explore interactions between DEGs. The GO analysis revealed that 1,992 DEGs were significantly enriched in 21 terms, including excretion, epidermis development, ion transmembrane transport, chemical synaptic transmission, chloride transmembrane transport, ion transport, potassium ion transmembrane transport, integral component of plasma membrane, extracellular region, extracellular space, plasma membrane, apical plasma membrane, anchored component of membrane, proteinaceous extracellular matrix, integral component of membrane, basolateral plasma membrane, calcium ion binding, heparin binding, sequence-specific DNA binding, transporter activity, and carbohydrate binding. In addition, the KEGG pathway analysis revealed that 1,992 DEGs were significantly enriched in five pathways, including neuroactive ligand-receptor interaction, calcium signaling pathway, gastric acid secretion, bile secretion, and pancreatic secretion. According to the STRING results, we constructed the PPI network. Then hub genes were selected with a high degree of interaction in the PPI network, including BDKRB2, NMUR2, NMU, BDKRB1, LPAR5, KNG1, LPAR3, SAA1, MCHR1, and PMCH. Further analysis of survival related to the expression of these hub genes revealed that BDKRB1, NMUR2, PMCH, and SAA1 are the key genes for the development of pRCC.
One hub gene, BDKRB1, is a well-established tumor suppressor gene, which is frequently mutated in familial breast and ovarian cancers. The gene product of BDKRB1 functions in a number of cellular pathways that maintain genomic stability, including DNA damage-induced cell cycle checkpoint activation, DNA damage repair, protein ubiquitination, chromatin remodeling, as well as transcriptional regulation and apoptosis. In this study, we found the role of BRCA1 in tumor suppression and DNA damage response, including DNA damage-induced cell cycle checkpoint activation and DNA damage repair. The other hub gene KNG1 (Kininogen-1) is expressed at low level in glioma cells. KNG1 can exert antiangiogenic properties and inhibit the proliferation of endothelial cells [11]. Previous work showed that KNG1 can be used as a serum biomarker for colorectal cancer [12]. Overexpression of the KNG1 inhibited proliferation and induces apoptosis of glioma cells [11]. In this study, KNG1 expression was downregulated in pRCC, which may be associated with the viability and angiogenesis of pRCC, but the analysis revealed no statistical impact of expression of this gene on survival, suggesting further investigation into the relationship between this gene and pRCC is required. Lysophosphatidic acid (LPA) is an extracellular biological lipid that interacts with G protein-coupled LPA receptors (LPAR1 to LPAR6) [13]. The lysophosphatidic acid receptor-3 (LPAR3) mediates viability among malignant cells and aggressiveness among certain tumors [14]. LPAR3 has been characterized as the major promoter of long-term viability in melanoma cells [15]. Other studies found that increased expression of LPAR3 increases malignancy in breast and ovarian cancers in vivo [16, 17]. In this study, LPAR3 was identified as a downregulated gene in pRCC. It was reported with the involvement of LPA5 in the activation of tumor progression in pancreatic cancer cells [13]. Bradykinin (BK) is produced in the inflammatory tissue microenvironment, where it acts in cell proliferation, leukocyte activation, cell migration, and endothelial cell activation [18]. BDKRB1 and BDKRB2 belong to the rhodopsin family of G protein-coupled receptors. The activation of BDKRB1 leads to the activation of macrophages, dendritic cells, and other cells in the tumor microenvironment, which have angiogenic properties and is related to the proliferation of malignant tumors [19]. BDKRB1 contributes to interleukin-8 production and glioblastoma migration [20]. Wang et al. reported that inhibition of BDKRB2, but not the B1 receptor, attenuated bradykinin-mediated invasion and migration in colorectal cancer cells and inhibited ERK1/2 activation and IL-6 production [21]. Thus, the identification of inhibitors against BDKRB1 may be a reasonable strategy to suppress pRCC. Neuromodulin U (NMU) activates the G protein-coupled receptor NMUR2, and NMU signaling interacts with several cancer-related pathways, including the WNT receptor cascade, resulting in increased activation of WNT/planar cell polarity (PCP) effector RAC1, which promotes tumor cell invasion and metastasis [22]. NMUR2 is a receptor that enhances NMU-mediated cell motility and invasion in human pancreas and endometrial cancer cells [23, 24]. Hub genes NMU and NMUR2 have not previously been reported to play roles in pRCC. PMCH encodes the 165 aa prohormone promelanin-concentrating hormone (PMCH), which is proteolytically processed into several peptides, including the oncogenic peptide melanin-concentrating hormone (MCH) [25]. In this study we found that increased expression of PMCH was associated with poor survival in patients with pRCC, suggesting PMCH may be a potential diagnostic biomarker or predictor of prognosis. Human serum amyloid A (SAA) is a high-density lipoprotein (HDL)-related lipoprotein with major roles in the regulation of inflammation and cholesterol transport [26]. Human serum amyloid A (SAA) has been widely regarded as an accurate and sensitive indicator of inflammation, which can be synthesized by the liver and cancer cells [27]. SAA1 regulates cell adhesion and migration and binding to laminin by inducing cytokine expression [28]. A previous study reported a relationship between increased SAA1 concentration and poor prognosis and distant metastasis in ccRCC patients [29].
In conclusion, bioinformatics analysis was used to identify DEGs that may be involved in the development or progression of the pRCC. This study identified several genes that may be involved in the pathology of papillary renal cell carcinoma. BDKRB1, NMUR2, PMCH, and SAA1 may contribute to the occurrence and development of papillary renal cell carcinoma. This identification of specific biological functions that may be involved in the mechanism of pRCC development provides new clues and directions for efforts to develop future treatments for papillary renal cell carcinoma.
Data Availability
All data generated or analyzed during this study are included within this article.
Conflicts of Interest
The authors declare that they have no conflicts of interests.
Authors’ Contributions
Le Li and XuKai Liu designed the manuscript. Pan Liu and Ting Sun transcribed the whole article. Yong Wen downloaded the data for analysis. All authors read and approved the final version of the submitted manuscript. Le Li, XuKai Liu, and Yong Wen contributed equally to this work.
Acknowledgments
This work was supported by Hunan Natural Science Foundation (2021JJ50068).