Key genes associated with prognosis and metastasis of clear cell renal cell carcinoma

Introduction

Clear cell renal cell carcinoma (ccRCC) is a highly aggressive malignancy that accounts for 85%–90% of all renal cell carcinoma cases (Chen et al., 2017). Approximately 17%–30% of ccRCC patients show distant metastasis at the time of diagnosis (Motzer, Bander & Nanus, 1996). Since ccRCC is not sensitive to radiotherapy and chemotherapy, radical or partial nephrectomy is still the main treatment, and 20% to 40% of cases that receive radical/partial nephrectomy still show recurrence or metastasis (Teng et al., 2014). Through the in-depth studies of ccRCC driving genes, studies have shown promising findings for molecular targeted treatments such as those targeting hypoxia-inducible factor 2α (HIF2α), indicating the potential success with new targeted therapies (Chen et al., 2016). However, treatments for patients with metastasis are still limited and the mechanism underlying ccRCC metastasis remains unclear. Hence, clarifying the mechanism on ccRCC metastasis is critical for of effective therapeutic targets to ccRCC patients with distant metastasis.

Microarray technology and bioinformatics analysis have facilitated the exploration of genetic alterations in various cancers. Bioinformatics analyses have helped researchers to identify genes specifically related with distant metastasis in a variety of tumors. For example, bioinformatics studies revealed that the expression of the sodium-potassium-chloride cotransporter 1 (NKCC1) gene was significantly higher in gastric cancer tissue than in normal gastric tissue. NKCC1 was shown to promote the proliferation, invasion, migration and epithelial to mesenchymal transition (EMT) of gastric cancer cells (Wang et al., 2021). Another study reported that U three protein 14a (UTP14A) stimulates the proliferation and metastasis of esophageal squamous cell carcinoma cells by regulating the PERK/eIF2a signaling pathway (Li et al., 2021).

With the aim to elucidate the key molecular events involved in metastatic ccRCC, we downloaded and analyzed three gene expression microarray datasets on ccRCC from GEO. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment of differentially expressed genes (DEGs) to investigate pathways involved in ccRCC metastasis. We also identified key factors associated with ccRCC prognosis and metastasis through protein-protein interaction (PPI) network construction, survival analysis and in vitro verification.

Materials and Methods

Results

Identification of DEGs in metastatic ccRCC

We searched for DEGs in metastatic ccRCC samples compared with primary ccRCC samples in three datasets. A total of 996 DEGs in GSE85258, 711 DEGs in GSE105288 and 20,820 DEGs in GSE2254 were obtained. Venn diagram analysis revealed 104 co-expressed DEGs among the three datasets (Fig. 1A).

PPI network construction, hub gene selection and key gene identification

Analysis using the STRING online database demonstrated the PPI network of the DEGs. The most significant module containing 17 hub genes was established by the MCODE app of Cytoscape (Figs. 2A, 2B). GO and KEGG pathway enrichment analyses of the 17 hub genes are shown in Table 1. GO enrichment results revealed that the 17 hub genes were mainly enriched in extracellular matrix organization, collagen fibril organization, collagen catabolic process, skeletal system development cellular response to amino acid, cell adhesion and protein heterotrimerization in BP; extracellular matrix, collagen trimer, endoplasmic reticulum lumen, proteinaceous extracellular matrix, extracellular region 14, extracellular space and extracellular exosome in CC; and extracellular matrix structural constituent, platelet-derived growth factor binding, integrin binding, heparin binding, SMAD binding 3, collagen binding and peptidase activator activity in MF. KEGG pathways were mainly enriched in protein digestion and absorption, ECM-receptor interaction, focal adhesion, PI3K-Akt signaling pathway, Amoebiasis and Platelet activation (Table 1).

Table 1:

GO and KEGG pathway enrichment analysis of 17 hub DEGs.

Term	Description	Gene count	Gene ratio (%)	FDR
GOTERM_BP_DIRECT
GO:0030198	extracellular matrix organization	12	70.59	2.39E−16
GO:0030199	collagen fibril organization	8	47.06	1.62E−13
GO:0030574	collagen catabolic process	8	47.06	4.33E−12
GO:0001501	skeletal system development	7	41.18	6.92E−08
GO:0071230	cellular response to amino acid stimulus	5	29.41	2.35E−06
GO:0007155	cell adhesion	8	47.06	2.35E−06
GO:0070208	protein heterotrimerization	3	17.65	0.00154
GO:0043588	skin development	3	17.65	0.00821
GO:0001568	blood vessel development	3	17.65	0.00912
GO:0009612	response to mechanical stimulus	3	17.65	0.01975
GO:1903225	negative regulation of endodermal cell differentiation	2	11.76	0.02424
GO:0048592	eye morphogenesis	2	11.76	0.03332
GO:0007229	integrin-mediated signaling pathway	3	17.65	0.04214
GOTERM_CC_DIRECT
GO:0031012	extracellular matrix	15	88.24	2.13E−22
GO:0005581	collagen trimer	9	52.94	5.36E−14
GO:0005788	endoplasmic reticulum lumen	10	58.82	1.33E−13
GO:0005578	proteinaceous extracellular matrix	10	58.82	2.06E−12
GO:0005576	extracellular region 14	14	82.35	4.64E−11
GO:0005615	extracellular space	12	70.59	4.99E−09
GO:0070062	extracellular exosome	10	58.82	7.93E−04
GO:1903561	extracellular vesicle	3	17.65	0.00302
GO:0005584	collagen type I trimer	2	11.76	0.00546
GO:0043202	lysosomal lumen	3	17.65	0.00614
GO:0005589	collagen type VI trimer	2	11.76	0.00614
GO:0005588	collagen type V trimer	2	11.76	0.00614
GO:0005796	Golgi lumen	3	17.65	0.00676
GOTERM_MF_DIRECT
GO:0005201	extracellular matrix structural constituent	7	41.18	7.00E−10
GO:0048407	platelet-derived growth factor binding	4	23.53	1.66E−06
GO:0005178	integrin binding	5	29.41	2.35E−05
GO:0008201	heparin binding	4	23.53	0.00310
GO:0046332	SMAD binding 3	3	17.65	0.00431
GO:0005518	collagen binding	3	17.65	0.00698
GO:0016504	peptidase activator activity	2	11.76	0.03911
KEGG_PATHWAY
hsa04974	Protein digestion and absorption	8	47.06	2.56E−11
hsa04512	ECM-receptor interaction	7	41.18	1.73E−09
hsa04510	Focal adhesion	7	41.18	2.12E−07
hsa04151	PI3K-Akt signaling pathway	7	41.18	3.37E−06
hsa05146	Amoebiasis	5	29.41	1.04E−05
hsa04611	Platelet activation	5	29.41	2.34E−05

DOI: 10.7717/peerj.12493/table-1

Notes:

GO

Gene Ontology

MF

Molecular Function

CC

cellular component

BP

biological process

KEGG

Kyoto Encyclopedia of Genes and Genomes

Four genes, PCOLCE, P4HB, COL6A2 and COL6A3 were sorted out with P-value < 0.01 among the 17 hub genes after survival analysis performed on the UALCAN website. Data showed that patients with high expression of PCOLCE, P4HB, COL6A2 and COL6A3 had worse survival than those with low expression of each gene (Figs. 2A, 2B). Also, these four genes were found to be higher expressed in primary ccRCC compared to normal kidney tissue (Fig. S1).

In vitro validation of the four key genes in metastatic ccRCC

qRT-PCR showed that PCOLCE, P4HB, COL6A2 and COL6A3 mRNAs were all upregulated in the Caki-1 cell line compared with its primary counterpart, the Caki-2 cell line (P < 0.05) (Fig. 3A).

We further examined the expressions of the proteins encoded by the four genes in 22 pairs of primary ccRCC and their corresponding metastatic counterpart tissue with immunohistochemistry. The TNM staging information of 12 cases of primary ccRCC at first diagnosis is presented in Table S1. Higher expressions of all four proteins were observed in most metastatic ccRCC tissues compared with the primary sites, as shown by mean IOD results (P < 0.05). Higher expressions of COL6A2 and PCOLCE in metastatic counterparts were observed in 19 out of 22 pairs, and higher expressions of COL6A3 and P4HB expression were detected in 17 out of 22 pairs (Fig. 3B).

Discussion

CcRCC is believed to originate from proximal tubular epithelium with aggressive behavior that was unparallel with its bland appearance. After surgery, the median survival time of ccRCC patients with metastasis is only 13 months, and the 5-year survival rate is 12.3% (Babaian et al., 2014). As one of the main genetic events contributing to ccRCC pathogenesis, the mutant VHL gene pathway has become a target for treatment towards ccRCC (Chen et al., 2016). With the ongoing effort in bioinformatic data management, novel candidate key genes may be identified as new therapeutic targets for ccRCC patients with metastasis. In order to obtain datasets with a long time span and sufficient samples from different institutions, we downloaded 3 datasets from GEO websites to screen out the intersection of differential expression genes.

In our study, we first obtained 104 co-expressed DEGs from three datasets. In the course of searching for key genes involved in ccRCC metastasis, we conducted bioinformatic analysis on 17 hub genes which were screened from the PPI network analysis to identify specific pathways that might be related to ccRCC metastasis. GO analysis showed the DEGs were mainly distributed in extracellular matrix, extracellular space, extracellular exosome and Golgi lumen of CC, and were mainly enriched in extracellular matrix organization, collagen fibril organization, cell adhesion and blood vessel development of BP, and extracellular matrix structural constituents, platelet-derived growth factor binding, integrin binding and collagen binding of MF. KEGG pathway analysis revealed that these DEGs were enriched in protein digestion and absorption, ECM-receptor interaction, focal adhesion and PI3K-Akt signaling pathway.

Exosomes, which functioned in carrying protein, ceramide, mRNA or micro-RNA as cargo, can induce metastasis through participation in EMT, mesenchymal-epithelial reverting transition and neoangiogenesis (Steinbichler et al., 2017). For example, in ccRCC cell lines, exosomes containing miR-19b-3p could initiate EMT and accelerate tumor metastasis by reducing PTEN expression (Wang et al., 2019a; Wang et al., 2019b). EMT was essential for tumor metastasis through inducing processes such as loss of cell polarity, disconnection with the basement membrane, acquiring plasticity, invasion and anti-apoptosis activities (Talbot, Bhattacharya & Kuo, 2012). Additionally, cell adhesion to ECM participates in multiple steps of tumor metastasis through ECM receptor interaction (Hamidi & Ivaska, 2018). For example, through the several binding sites from matricellular proteins members of the CCN family, integrin could help tumor cells resist anoikis and promote tumor cells metastasize (Lau, 2016; Simpson, Anyiwe & Schimmer, 2008). Studies have shown that the PI3K-Akt signaling pathway might induce the EMT process directly or through crosstalk with the Wnt/β-catenin signaling pathway, TGF-β or NF-kB and initiate distant metastasis (Xu, Yang & Lu, 2015). The phosphorylated AKT, a key kinase of the PI3K signaling pathway, was reported to be related to EMT in ccRCC (Wu et al., 2011). Our results also showed that DEGs in metastatic ccRCC were enriched in the PI3K signaling pathway, suggesting that PI3K signaling pathway–mediated EMT was likely to participate in the metastasis of ccRCC. Together, our results have revealed various processes that might be involved in ccRCC metastasis.

To explore the relationship between 17 hub genes and patient survival, the survival analyses were conducted. From survival analyses, PCOLCE, P4HB, COL6A2 and COL6A3 were found to be related to ccRCC patient survival, and these candidate genes were further tested with qRT-PCR and immunohistochemical staining in ccRCC cell lines and tissues. We found that the four genes showed higher expression in the metastatic counterpart of ccRCC both in cell lines and paired ccRCC tissues. P4HB and COL6A3 have been previously reported to be associated with ccRCC. And our study presents the linkage of PCOLCE and COL6A2 with ccRCC.

To the best of our knowledge, this is the first report that showed the relevance between PCOLCEand ccRCC. PCOLCE is a secreted glycoprotein that participates in the maturation of procollagen and ECM reconstruction by enhancing the activity of bone morphogenetic protein-1 (BMP-1) (Kessler, Mould & Hulmes, 1990; Moali et al., 2005; Zhu et al., 2019a; Zhu et al., 2019b). PCOLCE was overexpressed in osteosarcoma and promoted lung metastasis via twist family bHLH transcription factor 1 (TWIST1). Reducing PCOLCE expression could prevent the migration, invasion and lung metastasis of osteosarcoma cells (Wang et al., 2019a; Wang et al., 2019b). BMP-1 was reported to be up-regulated in gastric cancer and associated with poor survival and distant metastasis (Hsieh et al., 2018). Notably, a previous study showed that BMP-1 was elevated in ccRCC and promoted proliferation, migration and invasion in ccRCC cell lines (Xiao et al., 2020). Our study showed that the expression of PCOLCE, which functions as an enhancer of BMP-1, was higher in metastatic ccRCC. This finding raised the hypothesis that PCOLCE might be involved in ccRCC metastasis through regulation of BMP-1. Therefore, further research is needed to clarify the function and potential mechanism of PCOLCE in ccRCC metastasis.

P4HB, a beta subunit of prolyl 4-hydroxylase, helps cancer cells survive from apoptosis induced by endoplasmic reticulum stress (Lovat et al., 2008; Pajunen et al., 1991). P4HB promoted proliferation, invasion, migration and angiogenesis in glioma through the mitogen-activated protein kinase (MAPK) signaling pathway (Sun et al., 2017). In gastric cancer, P4HB was shown to play an important role in regulating invasion and migration in the HIF1α pathway (Zhang et al., 2019). Previous studies showed that HIF1α was activated in renal clear cell carcinoma and associated with a more aggressive phenotype (Maxwell et al., 1999; Na et al., 2003). And the overexpression of P4HB in ccRCC was associated with poor outcome (Zhu et al., 2019a; Zhu et al., 2019b). In our study, we found that P4HB was specifically elevated in metastatic ccRCC. Hence, these findings suggest that P4HB might play a key role in modulating ccRCC metastasis through HIF1α pathway.

In this study, we showed that COL6A2 and COL6A3 were related to ccRCC progression. COL6A1–3 are three subunits of collagen VI (Hessle & Engvall, 1984; Weil et al., 1988). Strong evidence has demonstrated the role of COL6A1–3 in promoting the progression of tumors (Chen, Cescon & Bonaldo, 2013). In serous ovarian cancer, COL6A2 was dysregulated and associated with poor prognosis through the TGFβ1 pathway (Cheon et al., 2014; Svoronos et al., 2020). COL6A3 functions in tumorigenesis and progression of cholangiocarcinoma through the E2F1/LMCD1-AS1/miR-345-5p/COL6A3 axis and serves as prognostic factor for pancreatic cancer (Yu et al., 2019). Up-regulation of COL6A1 in ccRCC was associated with poor prognosis of patients and promoted tumor growth in vivo (Wan et al., 2015). COL6A3 was reported to be overexpressed in metastatic renal cell carcinoma and associated with poor survival (Li et al., 2019). To date, no report has described the relation between COL6A2 and ccRCC. As COL6A2 and COL6A3 were identified as key genes in our results, further research should be carried out to examine their mechanisms in ccRCC metastasis.

This study has several limitations. First, we only examined 22 paired ccRCC samples and 2 cell lines to verify the bio-informatic findings. More samples need to be examined to confirm our results. Second, the mechanisms of the four key genes in regulating ccRCC cell metastasis are still unclear. Further exploration is required to clarify the function and potential mechanism of the four key proteins both in vitro and in vivo.

Conclusion

Our studies identified four key genes (PCOLCE, P4HB, COL6A2 and COL6A3) that were related to ccRCC patient survival and showed potential involvement in ccRCC metastasis. The proteins encoded by the four genes may be involved in the ECM-receptor interaction, focal adhesion and the PI3K-Akt signaling pathway. In vitro experiments confirmed increased expression of the candidate genes in metastatic ccRCC, both in cell lines and tissues. Future studies are required to evaluate the significance and mechanism of these key proteins in ccRCC.

Supplemental Information