Identification of Prognostic Markers of DNA Damage and Oxidative Stress in Diagnosing Papillary Renal Cell Carcinoma Based on High-Throughput Bioinformatics Screening

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.


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]. Tere 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]. Te molecular mechanism of pRCC has not been clearly defned. 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.
Te 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. Te 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 diferential gene enrichment (Gene Ontology, GO) analysis and KEGG pathway enrichment analysis, constructed PPI interaction networks, screened hub genes, and performed survival analysis.

Data Collection.
Te published transcriptome data related to papillary renal cell carcinoma were downloaded from TCGA (https://cancergenome.nih.gov/). Te data included 289 papillary renal cell carcinoma samples and 32 normal kidney tissues.

Identifcation 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 classifcation algorithms, including self-organized maps (SOM), artifcial neural networks (ANN), support vector machines (SVM), fuzzy logic, and hyphenated techniques such as neuro-fuzzy networks to standardize the data and analyze diferential expression. Genes with |logFC| > 2.0 and FDR <0.05 were considered diferentially expressed genes. To visualize the data graphically, the ggplot2 software package was used.

GO and KEGG Pathway Analysis.
Te DAVID database (DAVID; https://david.ncifcrf.gov) was used to perform annotation, visualization, and integrated discovery on the genes identifed as signifcantly diferently 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 P < 0.05 was considered statistically signifcant. Te 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 P value <0.05 was considered statistically signifcant.

Hub Genes Selection and Analysis of Modules from PPI
Networks. Te STRING database (http://string-db.org) aims to provide a critical assessment and integration of proteinprotein (PPI) interactions [8]. STRING was used to analyze the selected diferentially expressed genes and construct a PPI network. Ten, cytoHubba in Cytoscape software (version 3.7.2) was used to screen the top 10 hub genes in the PPI network [9].

Survival Analyses of Hub Genes.
Te expression profles and clinical data of 289 pRCC samples were downloaded from TCGA (http://tcga-data.nci.nih.gov) for the survival analysis of hub genes. Te Kaplan-Meier method was used for the survival analysis, and log-rank P values were calculated. A log-rank P value <0.05 was considered statistically signifcant.

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 confrmed 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.

Total RNA Isolation and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR). Te 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 quantifcations were normalized to the level of glyceraldehyde phosphate dehydrogenase (GAPDH) in the reaction. Te comparative threshold cycle (CT) method, which compares the diferences in CT values between common reference RNA and target gene RNA, was used to obtain the relative fold changes in gene expression. Te expressions were calculated by 2 −ΔΔct method. Each experiment was performed in triplicate and repeated three times.

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. Te 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 Chisquared 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. P < 0.05 was considered to be signifcantly diferent.

Identifcation of DEGs.
Te 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. Te data were normalized and logarithmized, probes without corresponding gene annotation information were removed, and repeated probes were removed to fnally get the expression profles 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 diferentially 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 signifcant diferences were selected (Figure 1(b)).

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. Te 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-specifc DNA binding, transporter activity, and carbohydrate binding. Te GO analysis fndings are shown in Figure 2 and Table 1.
We next performed KEGG pathway analysis to analyze the pathways at the functional level. Te results showed that DEGs were mainly enriched in neuroactive ligand-receptor interaction, calcium signaling pathway, gastric acid secretion, bile secretion, and pancreatic secretion. Te KEGG pathways associated with enriched DEGs associated with pRCC are presented in Figure 2(b) and Table 2.

Identifcation of Hub Genes and Analysis of Modules from PPI Networks.
Te 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

Survival Analysis of Hub Genes.
Expression data for a total of 289 pRCC samples were downloaded from TCGA. Te 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).

Te 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. Te results showed that the expression of hub genes BDKRB1, NMUR2, PMCH, and SAA1 in pRCC tissues was signifcantly 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 signifcantly 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 signifcant (Tables 3-6).    Figure 5). Ten we conducted the COX proportional risk model analysis. Te 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).

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. Te overall prognosis of pRCC is better than that of ccRCC, but pRCC prognosis is signifcantly worse than that of ccRCC when pRCC invades the renal vein and/or the inferior vena cava [10]. Tere is currently no specifc treatment for pRCC, and surgical treatment is the frst choice in clinical practice. Te prognosis of advanced patients is poor, a pRCC is insensitive to radiotherapy and chemotherapy. Terefore, 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   10 28 Journal of Oncology 7    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. Te 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    [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 proteincoupled LPA receptors (LPAR1 to LPAR6) [13]. Te 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 identifed 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 infammatory 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. Te 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]. Tus, the identifcation of inhibitors against BDKRB1 may be a reasonable strategy to suppress pRCC. Neuromodulin U (NMU) activates the G proteincoupled 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) efector 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 melaninconcentrating 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 infammation and cholesterol transport [26]. Human serum amyloid A (SAA) has been widely regarded as an accurate and sensitive indicator of infammation, 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. Tis study identifed 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. Tis identifcation of specifc biological functions that may be involved in the mechanism of pRCC development provides new clues and directions for eforts to develop future treatments for papillary renal cell carcinoma.

Data Availability
All data generated or analyzed during this study are included within this article.