Firstly, the expression level of NEIL3 in pan-cancer, including KIRC,was investigated by using the Tumor Immune Estimation Resource (TIMER.0) database (http://timer.cistrome.org/) (14). Level 3 RNA-sequencing data (HTSeq-Counts and HTSeq-FPKM) and clinical information of KIRC were downloaded from the UCSC Xena browser (https://xenabrowser.net/) (15). We converted HTSeq-FPKM data into transcripts per million (TPM) format for further analyses, and used HTSeq-Counts for differential analysis by the DESeq2 package. Cases with complete clinical characteristics and survival information were included in the study.

Microarray data of KIRC expression profiling (GSE40435, GSE53757 (16), GSE68417 (17) and GSE73731) were obtained from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (18). All gene expression microarray data were normalized using the limma R package. GSE40435, GSE53757 and GSE68417 were used to validate the differential expression of NEIL3 in clear cell carcinoma and normal tissues. GSE73731 with 265 KIRC samples was used to verify the relationship between NEIL3 expression and immune infiltration.

Clinical parameters of KIRC patients, including age, sex, histological grade, pathological stage, TNM stage, and vital status, were selected from previously downloaded clinical information for further analysis. Subsequently, based on median NEIL3 expression, we classified the patients into two groups. The survival curves were plotted to compare overall survival (OS), progression-free interval (PFI), and disease-specific survival (DSS) between two groups by the R package survminer (19). The receiver operating characteristic (ROC) curve was drawn by the R package pROC to verify the diagnostic value of NEIL3. We identified risk factors by univariate and multivariate Cox analysis. Furthermore, we construct a nomogram model by using the R package rms. The calibration curves were plotted to compare the consistency of the model prediction probability with the observed probability. The time-dependent ROC was performed by using the time ROC package to evaluate the prediction accuracy of the model. Moreover, we performed a decision curve analysis (DCA) to verify the clinical value of the predictive model.

Gene Ontology (GO) analysis was performed to predict NEIL3-related biological functions using the R package clusterProfiler (20). GSEA v4.2.3 software (21) was used for gene set enrichment analysis (GSEA) to identify NEIL3-related signal pathways. For each analysis, we set up 1000 permutations. The normalized enrichment score (NES) and false discovery rate (FDR) were applied to select significant pathway enrichment.

To further investigate the molecular mechanisms associated with NEIL3, we usedthe GEPIA browser (http://gepia.cancer-pku.cn/) (22) to identify co-expressed genes significantly associated with NEIL3. Subsequently, Spearman’s correlation analysis was carried out on the co-expressed genes. Meanwhile, the correlation between co-expressed genes and NEIL3 was also validated in the TIMER2.0 database. In addition, a protein-protein interaction (PPI) network was established through the STRING database (https://cn.string-db.org/) (23).

Gene Set Cancer Analysis (GSCA) (http://bioinfo.life.hust.edu.cn/GSCA/#/) is a comprehensive platform, which can provide methylation analysis (24). Thus, we used the data from GSCA online software to analyzed the relationship between NEIL3 methylation and survival in KIRC patients.

CIBERSORT was applied to calculate the proportion of 22 immune cells in KIRC patients (25). We performed single-sample Gene Set Enrichment Analysis (ssGSEA) to evaluate the relationship between NEIL3 expression and 28 immune cells by the R package GSVA (26). In addition, in order to predict the response of KIRC patients to immunotherapy, we also used TISIDB database (http://cis.hku.hk/TISIDB/) (27) to analyze the correlation between NEIL3 expression and immunomodulators.

All statistical analyses were performed by R (4.1.3). Chi-square or Fisher’s exact test was used to compare the clinical characteristics between high expression group and low expression group. Spearman’s coefficient was used for correlation analysis. The association between clinicopathological features and NEIL3 expression was analyzed using the Wilcoxon signed-rank sum test and logistic regression. In this study, all hypothetical tests were bilateral, and P < 0.05 was considered to be statistically significant.

Firstly, we assessed NEIL3 expression levels in different human tumors using the TIMER database. We found that NEIL3 was over-expressed in multiple tumors, including KIRC ( Figure 1A ). In addition, we further evaluated the expression levels of NEIL3 in KIRC using TCGA and GEO data, and found that the NEIL3 expression level in KIRC tissues was remarkably higher compared with that in normal kidney tissues (all P < 0.001) ( Figures 1B–E ).

To investigate the connection between clinical parameters and NEIL3 expression, our study incorporated clinical information of 539 KIRC samples from the TCGA database. Based on the median NEIL3 expression level, we divided these samples into two groups (low expression and high expression group). The results suggested that NEIL3 expression was closely related to histologic grade, pathologic stage, T stage, M stage, and vital status (all P < 0.001) ( Table 1 ). However, age (P = 0.698), gender (P = 0.953), and N stage (P = 0.062) had no correlation with NEIL3 expression ( Table 1 ).

Survival analysis showed that up-regulation of NEIL3 was relevant to worse OS, DSS, and PFI in KIRC patients ( Figures 2A–C ). Moreover, we performed ROC analysis to assess the diagnostic value of NEIL3 to identify KIRC tissues, and found that the area under curve (AUC) value was 0.944 ( Figure 2D ). The result indicated that NEIL3 had a strong diagnostic performance. In addition, to further explore the risk factors for patients with KIRC, we conducted univariate and multivariate analyses. Univariate Cox regression analysis revealed that age, histologic grade, pathologic stage, TNM stage, and NEIL3 expression were associated with OS ( Table 2 ). Multivariate Cox survival analysis showed that age (P = 0.046), histologic grade, distant metastasis (P < 0.001) and NEIL3 expression (P = 0.009) were independent predictors of poor prognosis in KIRC patients ( Table 2 ). Based on the multivariate cox regression analysis, a nomogram model was created to predict the probability of survival using NEIL3 expression levels, age, histologic grade, and M stage ( Figure 3A ). The calibration curves showed good agreement between the predicted OS and the observed OS at 1, 3, and 5 years ( Figure 3B ). The time ROC curve showed that the nomogram model could effectively predict the OS of KIRC patients at 1, 3, and 5 years ( Figure 3C ). Additionally, we obtained a consistency index (C-index) of 0.753. These results indicated that this model had a relatively reliable predictive performance.

In addition, we performed DCA analysis to verify the clinical practicability of the predictive model. The DCA curves revealed that the combined model had a higher net benefit in predicting OS at 1, 3, and 5 years compared to other models that included a single clinical variable ( Figures 3D–F ). The results indicated that the combination of NEIL3 with clinical features was beneficial to predict the prognosis of KIRC patients.

Next, we conducted a GO analysis to characterize the biological functions of NEIL3 in KIRC. The GO analysis showed that NEIL3 was engaged in multiple biological processes, including nuclear division, mitotic nuclear division, acute inflammatory response, kinetochore, chromosome condensation, serine hydrolase activity, and receptor-ligand activity signal pathway ( Figure 4A ). In addition, to further investigate NEIL3-related signal pathways in renal clear cell carcinogenesis, we performed GSEA analysis and found 16 signal pathways demonstrated significant differences (NES > 1.72, P-value < 0.05, FDR < 0.1) in the high expression group ( Table 3 ). Cell cycle, DNA replication, mismatch repair, and p53 signal pathway are the important signal pathways involved in tumorigenesis in the high NEIL3 group, while the signal pathways involved in the immune response included antigen processing and presentation, intestinal immune network for IgA production, primary immunodeficiency ( Figures 4B–H ).

In order to further explore the potential role of NEIL3 in renal clear cell carcinogenesis, we identified genes that were positively correlated with NEIL3 through the GEPIA database. The top ten related genes (r ≥ 0.73) are shown in the co-expression heatmap ( Figure 5A ). The interaction network for these genes was established by using the STRING database ( Figure 5B ). In addition, the correlation between NEIL3 and these genes was validated in the TIMER2.0 database. The analysis suggested that NEIL3 was significantly associated with KIFC1 (cor = 0.743, P = 1.66e-94), CDCA8 (cor = 0.79, P = 8.04e-115), TPX2 (cor = 0.786, P = 4.1e-113), ERCC6L (cor = 0.742, P = 2.57e-94), BIRC5 (cor = 0.726, P = 1.7e-88), KIF20A (cor = 0.731, P = 4.4e-90), GTSE1 (cor = 0.751, P = 6.11e-98), CENPA (cor = 0.781, P = 8.46e-111), ORC6 (cor = 0.653, P = 4.54e-66),AURKA (cor = 0.667, P = 7.06e-70) ( Figures 5C–L ).

Previous studies have illustrated that aberrant DNA methylation is closely related to the development and progression of multiple tumors (28, 29). In our study, we found that the methylation level of NEIL3 in renal clear carcinoma tissues was lower than that in normal renal tissues ( Figure 6A ). Spearman correlation analysis showed that there was a negative correlation between NEIL3 methylation and mRNA expression in KIRC (cor = −0.36, FDR = 6e-11) ( Figure 6B ). Furthermore, we found that compared with the hypermethylated patients, NEIL3 hypomethylated patients had worse OS, progression-free survival (PFS) and DSS, while there was not significantly difference in disease-free interval (DFI) between the two groups ( Figures 6C–F ). Overall, patients with NEIL3 hypomethylation (indicating higher expression of NEIL3 mRNA) had worse survival outcomes than patients with NEIL3 hypermethylation.

In order to explore the association between NEIL3 expression and immune infiltration, we performed CIBERSORT and ssGSEA. CIBERSORT analysis showed that the proportion of plasma cells, T cells CD8, follicular helper T cells (Tfhs), regulatory T cells (Tregs), and M0 macrophages was increased in the high expression group, whereas the fraction of resting NK cells, monocytes, and M2 macrophage was decreased in the high expression group ( Figure 7A ). Meanwhile, ssGSEA revealed that activated CD4 T cells, activated dendritic cells (DCs), activated CD8 T cells, macrophages, myeloid-derived suppressor cells (MDSCs), natural killer cells (NK cells), Tfhs and Tregs were highly expressed in the NEIL3 overexpression group ( Figure 7B ). Furthermore, the connection between NEIL3 co-expressed genes and immune infiltration in KIRC was verified by ssGSEA. It was found that the expression levels of NEIL3 co-expressed genes had a positive association with multiple immune infiltrating cells (including B cells, CD4 T cells, CD8 T cells, DCs, MDSCs, Tfhs, and Tregs, etc.) ( Supplementary Figure 1 ).

We further verified the relation between the expression of tumor-infiltrating lymphocytes (TILs) and NEIL3 via the TISIDB database. The results indicated that the abundance of activated CD4 T cells (r = 0.668, P < 2.2e-16), activated CD8 T cells (r = 0.422, P < 2.2e-16), DCs (r = 0.4, P < 2.2e-16), MDSCs (r = 0.339, P = 1.07e-15), memory B cell (r = 0.358, P = 8.74e-18), Tfhs (r = 0.332, P = 4.84e-15), NKT cells (r = 0.256, P = 2.2e-09), Type 1 T helper cell (Th1 cells) (r = 0.284, P = 3.04e-11), Type 2 T helper cell (Th2 cells) (r = 0.304, P = 8.55e-13) and Tregs (r = 0.231, P = 7.12e-08) was positively correlated with NEIL3 expression ( Figures 8A–J ). In addition, we divided 265 renal clear cell carcinoma samples from GSE73731 into two groups based on the median NEIL3 expression level. CIBERSORT and ssGSEA were used to assess the relationship between NEIL3 expression and the extent of immune infiltration. Similarly, we found a higher level of immune cell infiltration in the over expressed group ( Supplementary Figure 2 ).

Overall, the high-expression group of NEIL3 had stronger immune infiltration than the low-expression group, especially regarding CD4 T cells, CD8 T cells, Tfhs and Tregs.

Finally, we investigated the association between 24 immunoinhibitors and NEIL3 expression using the TISIDB database. The heatmap showed that more than half of 24 immunoinhibitors were positively related to NEIL3 expression in KIRC ( Figure 8K ). Among them, the expression of PD1, CTLA4, LAG3, TIGHT, IL10, and CD96 was significantly associated with NEIL3 expression (all spearman’s r ≥ 0.3, all P < 0.001) ( Figures 8L–Q ).