{'fraction_non_alphanumeric': 0.036086909755835214, 'fraction_numerical': 0.02538453264493923, 'mean_word_length': 5.124670619235837, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 2, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis as a novel prognostic biomarker associated with immune cell infiltration in KIRC

Long non-coding RNAs (lncRNAs) play a significant role in multiple human cancers as competing endogenous RNAs (ceRNAs). However, a systematic mRNA-microRNA (miRNA)-lncRNA network linked to kidney renal clear cell carcinoma (KIRC) prognosis has not been described. In this study, we aimed to identify the prognosis-related ceRNA regulatory network and analyzed its relationship with immune cell infiltration to predict KIRC patient survival. The MMP25-AS1/hsa-miR-10a-5p/SERPINE1 ceRNA network related to the prognosis of KIRC was obtained through bioinformatics analysis based on The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Meanwhile, we constructed a three-gene-based survival predictor model, which could be referential for future clinical research. Methylation analyses suggested that the abnormal upregulation of the SERPINE1 likely resulted from hypomethylation. Furthermore, the immune infiltration analysis showed that the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis could affect the changes in the tumor immune microenvironment and the development of KIRC by affecting the expression of chemokines (CCL4, CCL5, CXCL13, and XCL2). Tumor Immune Dysfunction and Exclusion (TIDE) analysis indicated that the high expression of SERPINE1 might be related to tumor immune evasion in KIRC. In summary, the current study constructing the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 ceRNA network might be a novel significant prognostic factor associated with the diagnosis and prognosis of KIRC.


INTRODUCTION
Renal cell carcinoma is a malignant tumor that originates from renal tubular epithelial cells, accounting for nearly 90% of renal malignancies and 3% of adult malignancies. 1 Kidney renal clear cell carcinoma (KIRC) is the most typical subtype of renal cell carcinoma. 2 It is usually asymptomatic, and approximately 30% of patients are diagnosed in the advanced stage. Advanced KIRC has an extremely poor prognosis due to its inherent resistance to radiotherapy and chemotherapy. 3 Currently, surgical resection is the most effective treatment for KIRC patients. 4 Therefore, a comprehensive understanding of the molecular mechanisms of KIRC and the development of effective early diagnostic and therapeutic strategies are urgently needed.
Competitive endogenous RNA (ceRNA) is a gene expression regulation mechanism proposed by Salmena et al. 5 in 2011. The ceRNA hypothesis holds that endogenous RNA molecules have microRNA (miRNA) target sites and can competitively bind to miRNAs, thereby indirectly regulating the expression of miRNA target genes. 6 This competitive miRNA binding effect is also called the miRNA sponge effect. A large number of studies have found that long non-coding RNAs (lncRNAs) can act as miRNA sponges to regulate target gene expression. 7,8 The lncRNA-miRNA-mRNA ceRNA network has been shown to play vital roles in progression and metastasis of multiple cancers, such as colorectal cancer, 9 pancreatic cancer, 10,11 liver cancer, as well as breast cancer. [12][13][14] Nevertheless, the pivotal lncRNA-miRNA-mRNA ceRNA networks that are significantly associated with prognosis of KIRC still need further research.
The tumor microenvironment (TME) is a mixture of fluids, immune cells, stromal cells, extracellular matrix molecules, and numerous cytokines and chemokines. 15 KIRC has the typical features of an immunogenic tumor, including numerous tumor-infiltrating T lymphocytes (TILs) and cytotoxic T cells, which identify and selectively destroy tumor cells, as well as circulating tumor-specific T cells. 16 However, KIRC can suppress the anti-tumor immunity of CD4 T cells, CD8 T cells, natural killer (NK) molecules, and dendritic cells through unclear mechanisms, successfully evading immune recognition. 17 The crosstalk between the tumor cells and tumor-infiltrating immune cells is usually modulated by the ceRNA networks. 18 Therefore, it is necessary to do better research on KIRC-infiltrating immune cells and ceRNA networks.
Advances in immunotherapies, particularly immune checkpoint blockade (ICB) and engineered T cells, have revolutionized cancer therapy in recent years. 19 ICB exerts profound anti-tumor effects in many cancer types. However, the efficacy of ICBs is greatly affected by the tumor microenvironment. 20 KIRC patients usually exhibit higher immune-related scores. Previous studies have shown that the interplay of somatic alterations and immune infiltration modu-lates the response of advanced KIRC to PD-1 blockade. 21 Although ICB has achieved durable disease control in some patients with advanced KIRC, the molecular mechanisms underlying ICB resistance have not been fully understood. 22 Therefore, it is necessary to do better research on the relationship between ICB response and ceRNA networks in KIRC.
In this study, we first obtained differentially expressed mRNAs (DEGs) by analyzing 4 GEO datasets and 1 The Cancer Genome Atlas (TCGA) dataset. We further identified and visualized 33 hub genes by a protein-protein interaction (PPI) analysis. Through prognostic analysis, 14 upregulated genes and 9 downregulated genes were used for subsequent analysis. The experimentally validated micro-RNA-target interactions database (miRTarBas) was used to predict the upstream miRNA of the above gene, and miRNet and starBase were used to predict the upstream lncRNAs of candidate miRNAs. Through differential expression analysis, co-expression correlation analysis, prognostic analysis, and Cox regression analysis, we successfully established a new ceRNA (MMP25-AS1/hsa-miR-10a-5p/SERPINE1) regulatory network ( Figure 1) and constructed a three-gene-based survival predictor model. Methylation analysis, immune infiltration analysis, and Tumor Immune Dysfunction and Exclusion (TIDE) analysis were further performed to study the potential biological function of the MMP25-AS1/hsa-miR-10a-5p/SER-PINE1 axis in KIRC. The established ceRNA network may help us to fully understand the pathogenesis of KIRC and provide KIRC with promising diagnostic biomarkers or effective therapeutic targets.

Identification of significant DEGs in KIRC
In the present study, five mRNA microarray datasets (GSE15641, GSE36895, GSE46699, GSE53757, and TCGA) were analyzed to obtain DEGs between KIRC tissues and adjacent normal tissues. As shown in the volcano plot, the genes with adjusted (adj.) p value < 0.05 and log2 (fold change) R 1.5 or % À1.5 were regarded as significant DEGs (Figure 2A). In the GSE15641 dataset, a total of 825 upregulated and 347 downregulated significant DEGs were screened out. In the GSE36895 dataset, there were 448 upregulated and 680 downregulated significant DEGs in KIRC tissues compared with adjacent non-tumor samples. For the GSE46699 dataset, there were 341 upregulated genes and 367 downregulated genes. For the GSE53757 dataset, a total 1,047 upregulated genes and 1,230 downregulated genes were identified. For the TCGA dataset, there were 345 upregulated mRNAs and 335 downregulated mRNAs. Venn diagram analysis revealed that there were 103 common upregulated DEGs and 58 common downregulated DEGs in five datasets ( Figures 2B  and 2C).
To explore the biological functions of the identified DEGs, we performed GO and KEGG pathway enrichment analysis using Metascape. GO analysis results indicated that for biological processes (BPs), the common significant DEGs were mainly enriched in response to hypoxia, response to wounding, and extracellular structure organization ( Figure S1A). For cell components (CCs), the  With regard to KEGG pathways, we found that the common significant DEGs were significantly enriched in PI3K-Akt signaling pathway, phagosome, and ECM-receptor interaction ( Figure S1D).
In order to understand the mutual interaction of the identified DEGs, PPI networks were constructed respectively for the upregulated and downregulated DEGs. According to the node degree calculated by CytoHubba of Cytoscape software, we obtained 33 hub genes, among which 20 were upregulated hub genes and 13 were downregulated hub genes ( Figures 2D  and 2E). Meanwhile, the overall survival (OS) of those hub genes in KIRC patients was also assessed by Kaplan-Meier (KM) plotter analysis with the TCGA database. We found that 14 upregulated hub genes (C1QA, C1QB, C3, CAV1,  CCL5, COL1A1, CXCR4, FN1, IL10RA, LOX,  MYC, SERPINE1, TYROBP, and VEGFA) were not only dramatically upregulated in KIRC but also significantly correlated with poor prognosis of KIRC patients ( Figure S2). On the other hand, 9 genes (ALB, ALDOB, AQP2, CLCNKB, G6PC, HMGCS2, PCK1, SCNN1A, and SLC12A1) were screened from 13 downregulated hub genes, and the low expression of them was related to the poor prognosis in KIRC ( Figure S3). The 23 key genes identified from the expression pattern and survival analysis were selected for next analyses.

Identification and validation of key miRNAs
The upstream miRNAs of those 23 hub genes were predicted by using the miRTarBase. Based on strong evidence validations (reporter assay, western blot, or qRT-PCR), a total of 112 upstream miRNAs (miRNA group1) were eventually identified as interacting with 10 upregulated hub mRNAs. The mRNA-miRNA regulatory networks for upregulated mRNA, consisting of 147 miRNA-mRNA relationships, was constructed and visualized by using the Cytoscape ( Figure 3A). Meanwhile, a total of 4 upstream miRNAs (miRNA group2) were eventually identified as interacting with 3 downregulated hub mRNAs. Cytoscape was used to construct and visualize a mRNA-miRNA regulatory network for downregulated mRNA, which consists of 6 miRNA-mRNA relationships ( Figure S4A). Furthermore, we screened the differentially expressed miRNAs between tumor and normal samples in KIRC. According to the cut-off criteria (adj. p value < 0.05 and log2 (fold change) R 1.5 or % À1.5), 76 miRNAs were differentially expressed between KIRC and normal tissues. As shown in the volcano plot, there were 15 upregulated miRNAs and 61 downregulated miRNAs ( Figure 3B). Venn diagram analysis showed that 17 of the 112 miRNAs targeted to the upregulated hub mRNAs were significantly downregulated in KIRC ( Figure 3C). However, among the 4 miRNAs that bound to the downregulated hub mRNAs, none of them were significantly upregulated in KIRC ( Figure S4B). Thus, we used the starBase database to analyze the expression correlation between 17 miRNAs (including 5 miRNAs and 5 mRNAs) and target mRNAs in KIRC (Table S1). We found 6 inversely correlated miRNA-mRNA pairs, including 5 miR-NAs and 5 mRNAs ( Figure 3D). Among the 5 miRNAs, 2 of them (hsa-miR-10a-5p and hsa-miR-200b-3p) were prognostic biomarkers according to the overall survival analysis for patients with KIRC (Figures 3E and 3F).

Construction of the lncRNA-miRNA-mRNA network
We first looked for those lncRNAs differentially expressed in KIRC patients in the TCGA database, generating volcano maps based on those genes with adj. p value < 0.05 and log2 (fold change) R 0.5 or % À0.5 ( Figure 4A). A total of 377 significantly differentially expressed lncRNAs were obtained, including 270 upregulated lncRNAs and 107 downregulated lncRNAs. Next, we predicted the upstream potential lncRNAs of the hsa-miR-10a-5p and hsa- miR-200b-3p through miRNet and starBase databases. The overlapping Venn diagram showed that 6 upregulated lncRNAs could potentially bind to hsa-miR-10a-5p, while 4 upregulated lncRNAs were predicted to bind to hsa-miR-200b-3p ( Figure 4B). Based on the above information, we constructed 10 miRNA-lncRNA interaction pairs ( Figure 4C). Subsequently, the expression correlation of miRNA-lncRNA interaction pairs in KIRC were evaluated by the starBase database (Table S2). A total of 2 significantly negatively correlated miRNA-lncRNA pairs were selected (Pearson's r < À0.1, p < 0.05) ( Figure 4D). At the same time, the two lncRNAs (MSC-AS1 and MMP25-AS1) showed a positive correlation with the corresponding mRNAs ( Figure 4E). Univariate Cox proportional hazards regression analysis indicated that MMP25-AS1, hsa-miR-10a-5p, hsa-miR-200b-3p, and SERPINE1 were significantly related to the overall survival of KIRC patients (p < 0.05) ( Figure 4F). Finally, we obtained a ceRNA network (MMP25-AS1/ hsa-miR-10a-5p/SERPINE1) related to the prognosis of KIRC, which is shown in Figure 4G. Furthermore, we explored the associations between gene expression (MMP25-AS1, miR-10a-5p, and SERPINE1) and clinicopathological characteristics, with the results summarized in Table 1. MMP25-AS1 expression was significantly associated with gender, pathological stage, American Joint Committee on Cancer (AJCC) T stage, and AJCC M stage. The expression of miR-10a-5p was significantly associated with race, pathological stage, AJCC T stage and AJCC M stage. SERPINE1 expression was significantly associated with gender, race, pathological stage, AJCC T stage, AJCC N stage, AJCC M stage, and tumor status. In order to further explore the biological function of the ceRNA network, KEGG pathway analysis of SERPINE1 via gene set enrichment analysis (GSEA) revealed that SERPINE1 was mainly involved in the chemokine signaling pathway, cytokine receptor interaction, ECM receptor interaction, and the JAK-STAT signaling pathway, which were closely related to epithelial-mesenchymal transition (EMT) ( Figure 4H).

Construction of a ceRNA-associated risk model
In order to explore KIRC's prognostic biomarker based on the ceRNA network (MMP25-AS1/hsa-miR-10a-5p/SERPINE1), we investigated the association of the ceRNA network and overall survival in KIRC patients using the univariate and multivariate Cox regression. KIRC patients were then divided into high-and low-risk subgroups according to the median cutoff of the three-gene-based prognosis risk score ( Figure 5A). Kaplan-Meier curves were plotted when the median risk score in the TCGA dataset was used as the cutoff value to compare survival risk between high-risk and low-risk groups. As shown in Figure 5B, KIRC patients with high risk scores had poor clinical outcomes. Receiver operating characteristic (ROC) analysis revealed that the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 prognostic model effectively predicted the survival of KIRC patients in 1-year, 3-year, and 5-year survival rates in the TCGA dataset ( Figure 5C).
Moreover, univariate analyses of clinical variables considered as potential predictors of survival are shown in Figure S5. Our results indicated that AJCC stage, T stage, age, recurrence, laterality, as well as risk score were significant risk factors for poor survival. Figure 5D shows the result of Cox multivariate analysis including these factors. However, multivariate analysis confirmed that AJCC stage, age, recurrence, and risk score were independent poor prognostic factor of overall survival. A nomogram was established based on the final Cox multivariable model, where the 1-, 3-, and 5-year survival probabilities of the KIRC patients could be easily predicted by summing up the points assigned to each respective variable ( Figure 5E).

Pan-cancer analysis of MMP25-AS1 expression and prognostic association
In order to examine the expression of MMP25-AS1 in cancers thoroughly, we conducted pan-cancer analysis of MMP25-AS1 expression in 23 types of cancers using the Gene Expression Profiling Interactive Analysis (GEPIA) database. The results showed that compared with normal tissues, the expression of MMP25-AS1 was significantly decreased in several types of tumor tissues, such as kidney chromophobe (KICH), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), thyroid carcinoma (THCA), thymoma (THYM), and uterine corpus endometrial carcinoma (UCEC). However, it was only significantly increased in KIRC ( Figure 6A). Later, univariate Cox analysis suggested that there were only significant associations between MMP25-AS1 expression and KIRC patient prognosis (Figure 6B). Meanwhile, analyses of the Clinical Proteomics Tumor Analysis Consortium (CPTAC) database revealed that the protein levels of SERPINE1 were significantly upregulated in KIRC tumor samples relative to normal tissues ( Figure 6C). The expression of MMP25-AS1 was associated with tumor stage, as it had a significantly higher expression level in KIRC stage III and IV compared to stage I patients ( Figure 6D). For this, all 510 patients were divided into an early-stage group (stages I and II) and an advanced-stage group (stage III and IV), and a survival analysis of each subgroup was conducted. We found that high MMP25-AS1 expression was significantly correlated with poorer survival in patients with early-stage and advanced-stage disease ( Figures 6E and 6F). Finally, we evaluated the effect of MMP25-AS1 expression on overall survival, progression-free survival (PFS), and disease-specific survival (DSS). MMP25-AS1 expression in KIRC tumor was significantly associated with overall survival , PFS, and DSS ( Figures 6G-6I).

Relationship between methylation and expression of SERPINE1
It has been shown that the occurrence and development of KIRC was associated with aberrant methylation of tumor suppressor genes. Our research found that DNA methyltransferases (DNMT1, DNMT3A, and DNMT3B), the key enzymes involved in DNA methylation, were significantly upregulated in the highrisk score group and the MMP25-AS1 high-expression group (Figure 7A). To further clarify the mechanism of abnormal upregulation of SERPINE1 in KIRC tissues, we explored the relationship between the expression level of SERPINE1 and its methylation status. The analysis of UALCAN demonstrated that the methylation level of SERPINE1 was significantly lower in KIRC tissues than in the normal tissues ( Figure 7B). Similarly, DiseaseMeth version 2.0 analysis by the Illumina 450K array showed a significantly decreased methylation level in KIRC tissues as compared with normal tissues ( Figure 7C). Additionally, we found the 3 most significant methylation sites (cg08792542, cg12584355, and cg01975495) in the DNA sequences of SERPINE1 that were negatively associated with their expression levels ( Figure 7D). Kaplan-Meier survival analysis based on methylation status showed that cg08792542 and cg01975495 site hypomethylation was associated with poor survival in KIRC ( Figure 7E).
Correlation between immune infiltration and the risk score of the ceRNA network In order to further study the correlation between the ceRNA network and tumor microenvironment immune cell infiltration, we first calculated the score of 28 immune cells' infiltration in each KIRC sample using ssGSEA. At the same time, the association between immune infiltration of each cell type and patient survival was assessed using Cox proportional hazards regression. A univariate Cox regression analysis showed that 15 immune infiltration cells were statistically significant for survival risk in KIRC cancer (Table S3). We used an unsupervised clustering method to divide all samples into three  Figure 8A). Survival analysis showed that patients in the high immune cell infiltration cluster had shorter overall survival than patients in the low immune cell infiltration cluster (Figure 8B). Further research found that the ceRNA network risk score was highly correlated to the immune cell infiltration. The ceRNA network risk score of patients in the high immune cell infiltration group was significantly higher than that of patients in the low immune cell infiltration group, and the expression of SERPINE1, hsa-miR-10a-5p, and MMP25-AS1 also showed the same trend (Figure 8C). Moreover, the ceRNA network risk score was negatively correlated with tumor purity ( Figure 8D). We found that ceRNA network risk score had significantly positive correlation with 10 cell

subsets (activated CD4+ T cells, activated CD8+ T cells, activated dendritic cells, macrophages, myeloid derived suppressive cells (MDSCs), effector memory CD4+ T cells, gamma delta T cells, natural killer T cells, central memory CD4+ T cells, and type 1 T helper cells)
, and the increase in these 10 cell subsets was associated with poor prognosis (Figures 8D and 8E). Only CD56bright natural killer cells had significant negative correlation with the ceRNA network risk score, and lower CD56bright natural killer cells had a poor prognosis ( Figures 8D and 8E).
The high expression of SERPINE1 might be related to tumor immune evasion Immune checkpoints were negative regulators of T cell activation, T cell proliferation, and effector functions. ICB was a promising approach to activating antitumor immunity. Hence, we evaluated 8 genes previously reported to be targets of immune checkpoint inhibitors: CD274 (PD-L1), CTLA4, HAVCR2, LAG3, PDCD1 (PD-1), PDCD1LG2 (PD-L2), TIGIT, and SIGLEC15. We found that 7 genes other than SIGLEC15 were significantly highly expressed in KIRC and were significantly related to tumor stage ( Figure 10A). The relationship between SERPINE1 expression and the expression of immune checkpoint genes was analyzed using the Spearman correlation analysis. The results showed that SERPINE1 expression level was significantly positively correlated with 4 immune markers (CTLA4, LAG3, PDCD1LG2, and TIGIT) (|cor| > 0.1 and p < 0.01) ( Figure 10B).
Afterward, the TIDE algorithm was used to predict the likelihood of an immunotherapy response. A low TIDE score meant a high probability of response to ICB therapy. The result showed that KIRC patients with early-stage disease (AJCC stage I and II) had lower TIDE scores compared to the patients with advanced-stage disease (AJCC stage III and IV) ( Figure 10C). Meanwhile, we observed that TIDE scores were significantly lower in the SERPINE1 low-expression group compared with the SERPINE1 high-expression group (Figure 10D). The above results indicated that KIRC patients with high expression of SERPINE1, regulated by the MMP25-AS1/hsa-miR-10a-5p axis, might have a higher chance of anti-tumor-immune escape and thus showed a lower ICB treatment response rate.

DISCUSSION
KIRC is a common urinary system tumor with a high level of tumorinfiltrating immune cells and high aggressiveness. Because KIRC is insensitive to conventional radiotherapy and chemotherapies, treatment of this disease largely relies on targeted therapy and immunotherapy. Due to the lack of precise and specific targets, the longterm efficacy of molecular targeted therapy and immunotherapy in the treatment of advanced KIRC is still not ideal. 23 The ceRNA regulatory network has been demonstrated to be involved in the initiation and progression of human cancers. 24 To the best of our knowledge, a comprehensive KIRC ceRNA network based on immune cell infiltration and prognostic factors has not been constructed so far. Therefore, we tried to establish the mRNA-miRNA-lncRNA ceRNA network related to the immune infiltration and prognosis of www.moleculartherapy.org KIRC and provide clues for further exploring KIRC prognostic biomarkers and potential therapeutic targets.
Here, for the first time, we searched for upstream lncRNAs and miR-NAs based on the key mRNAs and finally obtained the ceRNA network (MMP25-AS1/hsa-miR-10a-5p/SERPINE1) related to the immune infiltration and prognosis of KIRC. By searching for these genes in PubMed, we observed that miR-10a-5p and SERPINE1 had been studied for their roles in cancer or in the relationship with KIRC cancer. hsa-miR-10a-5p, as an antitumor gene, was significantly downregulated in KIRC specimens. Simultaneously low expression of hsa-miR-10a-5p was associated with poor prognosis. 25,26 SERPINE1 encoded plasminogen activator inhibitor-1 (PAI-1), which could modulate the essential processes of tumor development, growth, invasion, and metastasis, as well as angiogenesis and fibrosis. [27][28][29] We found that SERPINE1 (also known as PAI-1) was significantly upregulated in KIRC, and high expression of SERPINE1 was associated with poor prognosis. These results were consistent with previous studies (Sui et al., 2021 30 ). In this study, we first identified the MMP25 antisense RNA 1 (MMP25-AS1) as an oncogenic factor in KIRC. MMP25-AS1 expression was significantly higher in KIRC than in normal tissues, and high expression of MMP25-AS1 was associated with poor outcome by survival analysis. However, so far no studies have investigated the role of MMP25-AS1 in cancer.
DNA methylation is an epigenetic modification that maintains gene silencing with the addition of methyl groups to the fifth carbon position of the cytosine residues by DNA methyltransferases (including DNMT1, DNMT3A, and DNMT3B). 31 Aberrant DNA methylation, in particular promoter hypermethylation and transcriptional silencing of tumor suppressor genes, has an important role in the development of many human cancers, including KIRC. 32 In the present study, we demonstrated for the first time that DNA methyltransferases were highly expressed in KIRC patients with high expression of MMP25-AS1 and high ceRNA network risk scores. This suggests that MMP25-AS1 may be related to cancer-related DNA hypermethylation in KIRC patients. Meanwhile, we found that the level of SERPINE1 methylation was significantly lower and the expression of SERPINE1 was significantly higher in KIRC tissues compared to normal tissues, suggesting that the expression of SERPINE1 was upregulated by low DNA methylation in KIRC. Upon further research, we found that high expression of SERPINE1 was negatively correlated with methylation sites cg08792542 and cg01975495, and these methylation sites were associated with better survival in KIRC.
Immune cells that infiltrate tumors form an ecosystem in the tumor microenvironment to regulate cancer progression and are closely associated with clinical outcome in KIRC. 17 It has been found that KIRC tissue is often infiltrated by lymphocytes (T cells, B cells, and natural killer cells), macrophages, and dendritic cells. 33,34 The activated CD8+ T cells and activated dendritic cells were associated with the poor prognosis of KIRC, and the infiltrating CD4+ T cells could enhance KIRC cell proliferation. Infiltrating macrophages increased KIRC cell invasion via alteration of EMT and increasing the stem-cell-like population. 33 Thus, KIRC is considered to have a unique TME, because CD8+ T cell infiltrates and a high density of activated dendritic cells correlate with favorable prognosis in the majority of solid tumors except for KIRC. 35 In this study, we first studied the correlation between the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis risk score and 28 immune cells' infiltration. We found 11 types of tumor-infiltrating immune cells (activated CD4+ T cells, activated CD8+ T cells, activated dendritic cells, macrophages, MDSCs, effector memory CD4+ T cells, gamma delta T cells, natural killer T cells, central memory CD4+ T cells, type 1 T helper cells, and CD56bright natural killer cells) were not only highly correlated with the ceRNA network risk score but also particularly relevant to the prognosis of KIRC patients. These findings together suggest that these differences induced by the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis may have an impact on the changes in the tumor immune microenvironment and the development of KIRC.
Chemokines are chemoattractant cytokines that play a pivotal role in regulating migration and infiltration of immune cell populations. 36 Over 40 human chemokines are now acknowledged, each with its own specific pattern of cellular chemotaxis. 37 We found that only 4 of the 40 chemokines (CCL4, CCL5, XCL2, and CXCL13) were significantly correlated with the expression of the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis and 7 immune infiltrating cells (activated CD4+ T cells, activated CD8+ T cells, activated dendritic cells, macrophages, MDSCs, central memory CD4+ T cells, and type 1 T helper cells). Previous studies have reported that the RNA transcript levels of CCL4, CCL5, and CXCL13 are higher in KIRC tissues compared with normal tissues. 38,39 CCL5 is expressed not only in immune cells but also in tumor cells. Several studies have focused on the effect of CCL5 on tumors, finding that CCL5 can significantly promote tumor growth, metastasis, 40 angiogenesis, 41,42 and immune escape. 43,44 The experimental results show that reduced expression of CCL5 and CXCL13 decreased cell proliferation and invasion in the KIRC cell lines. 45,46 In addition, the latest study found that CXCL13 secreted by CD8+ T cells impairs the immune function of the infiltrating CD8+ T cells in KIRC. 47 At the same time, CCL5 and CXCL13 have been shown to be potential biomarkers and therapeutic targets, which is related to CD8+ T cell infiltration in KIRC. 46,48 The chemokine XCL2 was recently shown to play a critical role in recruiting cross-presenting dendritic cells to tumors. 49 In addition, studies have reported that tumor-resident natural killer cells express many chemokine genes (such as CCL4, CCL5, and XCL2), which are important for the infiltration of dendritic cells, T cells, and other immune cells. 50,51 These results indicated that the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis may participate in immune cell infiltration through chemokines CCL4, CCL5, XCL2, and CXCL13.
Immune checkpoints play important roles in immune regulation, and blocking immune checkpoints on the cell membrane is a promising strategy in the treatment of cancer. 52 Monoclonal antibodies targeting a variety of immune checkpoint inhibitors (e.g., anti-CTLA4, anti-PDCD1, and anti-LAG3) have been investigated across multiple tumor types, including prostate cancer, lung cancer, and KIRC. 22 However, many advanced KIRC patients are resistant to ICB therapy. 53 Through TIDE analysis, we found that the high expression of SERPINE1 in KIRC patients may be related to the resistance to ICB therapy. KIRC patients with high expression of SERPINE1 had high expression of multiple immune checkpoints (CTLA4, LAG3, PDCD1LG2, and TIGIT). Therefore, single-agent immunotherapy is usually inadequate, and combined immunotherapy may be a good choice for the treatment of KIRC cancer. Some studies have shown that blockade of PD-1 combined with blockade of CTLA4, LAG3, or TIGIT, to some extent, can reverse T cell dysfunction and enhance antitumor immunity. 54,55 The anti-CTLA4 antibody ipilimumab, in combination with the anti-PD-1 antibody nivolumab, has been US Food and Drug Administration (FDA) approved for the treatment of metastatic melanoma, metastatic colorectal cancer, and advanced KIRC. [56][57][58] PD-L1 is encoded by CD274 and PD-L2 is encoded by PDCD1LG2, both of which are immunosuppressant proteins inhibiting cytokine production and cytolytic activities of CD4+ and CD8+ T cells. 59,60 Our study also found that PD-L2 expression was significantly positively correlated with SERPINE1 expression but not PD-L1 expression. These results could further explain the mechanism of CD8+ T cell infiltration and poor prognosis of KIRC in renal cancer.
By integrating these mRNA-miRNA and miRNA-lncRNA interactions, a potential mRNA-miRNA-lncRNA ceRNA network associated with prognosis and immune cell infiltration of KIRC was constructed. However, there are some limitations to this study. First, the function of the key lncRNA MMP25-AS1 in KIRC requires further experimental research. Second, the ceRNAs axis of MMP25-AS1/hsa-miR-10a-5p/SERPINE1 should also be verified in in vitro and in vivo studies.

Conclusion
In conclusion, through integrated identification and means of stepwise reverse prediction, we identified KIRC-related lncRNA, miRNA, and mRNA. Based on the rules of ceRNA hypothesis, a ceRNA regulatory network consisting of the MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis was constructed successfully, and each component in the network was significantly associated with the prognosis and immune cell infiltration of KIRC patients. The MMP25-AS1/hsa-miR-10a-5p/SERPINE1 axis identified in our research provided valuable clues for further basic and clinical research and provided guidance for future KIRC diagnosis, targeted therapy, and immunotherapy.

Screening of differentially expressed genes
We identified DEGs, miRNAs (DEmiRNAs), and lncRNAs (DElncR-NAs) using the Linear Models for Microarray data (Limma) package 61 in R. The DEGs and DEmiRNAs were identified using the adj. p value < 0.05 and log2 (fold change) R 1.5 or % À1.5. We determined the DElncRNAs with thresholds of |log2 (fold change)| R 0.5 and adj. p value < 0.05. Meanwhile, volcano plots of the DEGs, DEmi-RNAs, and DElncRNAs were visualized using the R gplots program.
To investigate the potential function of the DEGs in the KIRC, GO functional enrichment analyses of the identified DEGs were performed using Metascape. 62 The identified DEGs were classified in terms of the BP, MF, and CC categories.

Identification of hub genes
The PPI network was constructed separately for upregulated and downregulated DEGs using the search tool for the retrieval of interacting genes online database 63 (https://string-db.org/). The PPI pairs with a combined confidence score R 0.4 were visualized in the network. The hub genes in the PPI network that had high degrees of connectivity were identified using the CytoHubba plugin for Cytoscape software (version 3.8.0). Finally, the top 20 upregulated DEGs according to the degree algorithm were selected separately, and the top 13 downregulated DEGs were selected separately as hub genes for further analysis. The upstream miRNAs that interacted with hub mRNAs were predicted by using the miRTarBase database 64 (https://mirtarbase. cuhk.edu.cn/). Targets of miRNAs were downloaded from the miRNA targets database miRTarBase, and only the miRNA/target pairs with strong experimental evidence were retained. We used miRNet database (https://www.mirnet.ca/) and starBase database 65 (http://starbase.sysu.edu.cn/index.php) to predict the upstream lncRNAs of miRNA. The common genes were selected for further analysis.

Co-expression analysis
The co-expression analysis of miRNA-mRNA, miRNA-lncRNA, and mRNA-lncRNA pairs identified in this study was assessed in KIRC cohorts from the starBase database. Gene pairs with |r| > 0.1 and p value < 0.05 were considered as potential pairs and were selected for further analysis.
lncRNA-miRNA-mRNA network construction Based on the negative correlation between miRNA and its targets, and the positive correlation between mRNA and the lncRNA/pseudogene, we constructed a lncRNA/pseudogene-miRNA-mRNA regulatory network. To further explore the molecular mechanism underlying the pathogenesis of KIRC, KEGG analysis of mRNAs in the network was performed by using the GSEA software package. Only the terms with a nominal p value < 0.05 and false discovery rate (FDR) q value < 0.25 were considered statistically significant.

Survival analysis and construction of a specific prognosis model for KIRC
A standard Cox proportional hazards model implemented in the R package survival was used for patient survival and Kaplan-Meier plotting. The time-dependent ROC curves were constructed using the survivalROC package of R. We developed a prognostic factor-based risk stratification nomogram for 5-year overall survival with Cox proportional hazards regression analysis using the rms library in R.

Immune infiltrate levels and expression analysis of ceRNA network
The marker genes of 28 immune cell types for ssGSEA were obtained from a previous study (Bindea et al., 2013 69 ). Infiltration levels for different immune cell types were quantified using the ssGSEA implementation R package gsva. Tumor purity, which represents the heterogeneity of each tumor sample, was estimated by the ESTIMATE R package. Heatmaps of immune infiltration levels were generated using the R heatmap package. Visualization scatterplot for a Spearman correlation between the ceRNA network risk scores and tumor-infiltrating immune cells in KIRC could be generated through the ggplot function in the ggplot2 package in R software. Correlation plots of all variables were produced using the corrplot function of the corrplot R package.

TIDE analysis
The TIDE model was a computational method that integrated the expression signatures of T cell dysfunction and T cell exclusion to model tumor immune evasion. The TIDE score of patients with KIRC from the TCGA dataset were downloaded from the TIDE website (http://tide.dfci.harvard.edu) after uploading the transcriptome profiles.

Statistical analysis
Most of the statistical analysis has been done by the bioinformatic tools mentioned above. The R software (version 4.0.4) was used for all the rest of the statistical analyses. Differential expression levels of mRNAs, miRNAs, and lncRNAs were estimated by a two-tailed Student's t test. The Benjamini-Hochberg FDR method was used for p value adjustment. Fisher's test was used to identify the significant GO terms. Spearman correlation coefficients were calculated to evaluate the correlations. A p value < 0.05 was considered as statistically significant.

DECLARATION OF INTERESTS
The authors declare no competing interests.