The mRNA expression profile data and corresponding clinical information of 539 KIRC samples and 72 adjacent non-tumor samples were acquired from TCGA website, and the 333 SARS-CoV-2 infection–related proteins (Supplementary Table S1) were summarized in the Human Protein Atlas (HPA) database. Then, we utilized the “limma” R package to identify the differentially expressed genes (DEGs) between KIRC tissues and adjacent normal tissues in the TCGA-KIRC cohort and calculated the intersections between DEGs and SARS-CoV-2–related genes using an online tool Venny2.1. The threshold for significance was set at a false discovery rate (FDR) < 0.05 and |log2fold change (FC)|  ≥ 1. Besides, the interaction network for overlapping SARS-CoV-2–related DEGs was generated by the “GeneMANIA” online database to predict their function (Franz et al., 2018).

Based on SARS-CoV-2–related DEGs, the “survival” R package was used to construct a prognostic model. The KIRC patients in the SARS-CoV-2–related gene signature building and validation need to meet the following inclusion criteria: 1) complete SARS-CoV-2–related DEG expression profile; 2) complete clinical data (including age at diagnosis, gender, TNM stage, histological grade, survival status, and survival time); 3) the follow-up or survival time must be longer than 30 days. After that, 501 KIRC patients were included for the subsequent signature construction in total. We randomly divided the TCGA-KIRC dataset (n = 501) into a training cohort (n = 252) and a testing cohort (n = 249) to construct and validate the prognostic signature, respectively.

Next, we performed a univariate Cox regression analysis to evaluate the prognostic value of SARS-CoV-2–related DEGs in the training cohort and screen genes significantly associated with the overall survival (p < 0.01). Subsequently, the multiple stepwise Cox regression analysis was applied to develop an optimal predictive model for the overall survival on the preliminary-screened prognosis-related genes. According to the smallest Akaike information criterion (AIC) value, the regression coefficients of each gene in this signature were generated (Song et al., 2020). The risk score for each KIRC patient was calculated as follows: Risk score= ∑ i n (Expression of  gene i ×Coefficient i) . The patients were ordered according to their risk scores and categorized into the high-risk and low-risk groups based on the median risk score.

The Kaplan–Meier method was used to plot survival curves and compare the differences of survival between high- and low-risk groups in the training cohort (n = 252). Additionally, the time-dependent receiver operating characteristic (ROC) curve analyses were implemented to assess the predictive capability for 1-, 3-, and 5-year OS of the prognostic signature and other clinical parameters by using the “survivalROC” R package. We also performed univariate and multivariate analyses with Cox regression to estimate the effect of the gene signature together with other clinical characteristics on survival (Long et al., 2019).

Meanwhile, the robustness and accuracy of the gene signature were further validated in the testing cohort and the entire TCGA-KIRC cohort. Besides, the E-MTAB-1980 dataset included 101 KIRC patients with reliable prognostic information which was downloaded from the ArrayExpress database for external validation (Sato et al., 2013).

The online database GEPIA was utilized to confirm the differential mRNA expression of the signature constituent SARS-CoV-2–related genes based on GTEx and TCGA datasets, and the protein expression was detected in the Human Protein Atlas (HPA) database (Tang et al., 2019; Huang et al., 2020). The association of the signature constituent genes with OS was further verified in the Kaplan–Meier Plotter website, and the genetic alteration analyses were performed through the cBioportal database (Gao et al., 2013; Nagy et al., 2021). Besides, the relationships between the mRNA expression levels of risk signature genes and six immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and Dendritic cells) infiltration were investigated by using the TIMER database (Li et al., 2017).

PCA scatter plots were generated to profile the distribution patterns of KIRC patients in two groups via using the “scatterplot3d” R package. Gene set enrichment analysis (GSEA) was achieved by GSEA v4.0.3 software to detect different functional phenotypes between the high- and low-risk groups based on GO Biological Process ontology gene sets. FDR <0.25 and a nominal p-value < 0.05 were considered to be significant.

We used the “CIBERSORT” R package to analyze the profiles of 22 tumor-infiltrating immune cells and quantify the relative proportions in the entire KIRC cohort (p < 0.05). Then, we further analyzed the correlations between the 22 types of immune cells (Newman et al., 2015). The “ESTIMATE” algorithm was also applied to calculate the immune and stromal scores and infer the tumor purity of each KIRC sample using gene expression data (Yoshihara et al., 2013). In addition, the correlations between the risk score signature and five immune checkpoints (CTLA-4, LAG-3, TIGIT, TIM-3, and PDCD1) were determined by the Pearson correlation test, and their expression levels were compared between high- and low-risk patients.

The clinical parameters including age, gender, TNM stage, histological grade, and the risk score signature were included in the development of the predictive nomogram to predict the overall survival of KIRC patients in the training test via the “rms” R package (Balachandran et al., 2015). The training test was split into two groups: high- and low-score groups, based on the total points of the nomogram. Then, we generated the Kaplan–Meier survival curve, ROC curve, and calibration plots of 1,3,5-OS to assess the predictive performance of this prognostic nomogram. Additionally, the good prediction performance of the nomogram was also further validated in the testing set, entire KIRC cohort, and E-MTAB-1980 dataset.

All statistical analyses were conducted by R software 3.6.3. The survival comparisons were carried out through the Kaplan–Meier method with the log-rank test. The p < 0.05 was considered to be statistically significant, unless otherwise stated.

The TCGA-KIRC dataset contains RNA sequencing data of 539 tumors and 72 adjacent normal tissues, and a total of 3,626 DEGs (2,513 upregulated and 1,113 downregulated) were identified between KIRC and normal kidney tissues (Figure 2A). We obtained the SARS-CoV-2–related gene list which contains 333 genes from the HPA online database and acquired 31 overlapped genes (SARS-CoV-2 related DEGs) including 15 upregulated and 16 downregulated between all DEGs and SARS-CoV-2–related genes by applying the Venn diagram (Figure 2B). The volcano plot and heatmap showing 31 SARS-CoV-2–related DEGs in Figures 2C,D. Moreover, the protein–protein interaction (PPI) network of SARS-CoV-2–related DEGs constructed using GeneMANIA uncovered these genes were primarily responsible for extracellular matrix and structure organization, exosome, and exonucleolytic nuclear-transcribed mRNA catabolic process, and deadenylation-dependent decay (Supplementary Figure S1).

In total, 501 KIRC patients were eligible for our inclusion criteria and included in the risk score model construction and validation. All of them, 252 patients were randomized to the training set for model construction, and 249 were randomized to the training set for model validation. The detailed clinical information of these patients is summarized in Table 1. Through univariate Cox regression analysis, we detected 11 SARS-CoV-2–related DEGs significantly associated with prognosis (p < 0.01) in the training set (Supplementary Figure S2A). Subsequently, the stepwise multivariate Cox regression analysis was performed based on above 11 prognostic SARS-CoV-2–related genes, and a 5-gene prognostic signature was constructed in the training set (Supplementary Figure S2B; Supplementary Table S2). The detailed description about Covid-19 bait, tissue specificity, blood specificity, and subcellular location of the five hub SARS-CoV-2–related genes (ACADM, CENPF, KDELC1, PLOD2, and TRMT1) in this predictive model is listed in Supplementary Table S3. Among these five genes, only downregulated ACADM was the protective factor with Cox HR < 1, and the other upregulated CENPF, KDELC1, PLOD2, and TRMT1 were unfavorable factors with Cox HR > 1. According to the expression and multivariate Cox regression coefficients of these five genes, we can calculate the risk score of each KIRC patient as follows: Risk score = (−0.4955 * ACADM Exp) + (0.3805 * CENPF Exp) + (0.2904 * KDELC1 Exp) + (0.2048 * PLOD2 Exp) + (0.4925 * TRMT1 Exp).

Then, the 252 KIRC patients in the training set were separated into low-risk (n = 126) or high-risk (n = 126) subgroups, using the median of the 5-gene signature risk scores as the cut-off (Figure 3A). The scatter plot for survival time and survival status distribution of each patient in the training cohort is illustrated in Figure 3B; patients with high risk usually had a higher likelihood of death earlier than those with low risk. Furthermore, the expression profile of the five genes in this signature was displayed as a heat map (Figure 3C). The Kaplan–Meier survival curves of two groups demonstrated that high-risk patients had a significant shorter survival time than those patients in the low-risk group, with p < 0.001 (Figure 3D). Moreover, the area under the curve (AUC) values of time-dependent ROC curves at 1-, 3-, and 5-years reached 0.766, 0.705, and 0.770, respectively; the results indicated that our prognostic signature had an excellent predictive ability (Figure 3E). Besides, the multivariate Cox regression analyses of the risk score signature and other clinical characteristics (age, gender, grade, and stage) were performed in the training set, and the risk score signature (HR = 1.302, 95% CI = 1.196–1.417, p < 0.001) was considered as an independent prognostic indicator for OS (Figures 3F,G; Supplementary Table S4).

To further examine the predictive ability of this risk score model, we also calculated the risk scores of the KIRC patients in the testing cohort (n = 249) using the same formula (Figure 4). The Kaplan–Meier survival analysis for two groups demonstrated that high-risk patients had a significant poor prognosis than patients with low risk, with p < 0.001 (Figure 4D). And, the AUC values of time-dependent ROC curves at 1-, 3-, and 5-years reached 0.741, 0.715, and 0.712, respectively (Figure 4E). Besides, the results of multivariate Cox regression analyses revealed that the risk score signature (HR = 1.142, 95% CI = 1.072–1.218, p < 0.001) was an independent prognostic indicator for OS in the testing cohort (Figures 4F,G; Supplementary Table S4). In addition, we also obtained the similar results in the entire TCGA-KIRC cohort (n = 501) (Figure 5).

The mRNA and protein expression levels of the five hub genes (ACADM, CENPF, KDELC1, PLOD2, and TRMT1) were validated in GEPIA and HPA databases. The mRNA expression level of ACADM was significantly reduced in KIRC tissues, while those of CENPF, KDELC1, PLOD2, and TRMT1 were significantly increased in KIRC tissues, which was determined by GEPIA (Figure 6A). The protein expression levels were detected through immunohistochemistry analyses in the HPA database, and the immunohistochemistry results further confirmed their expression difference between normal kidney tissues and KIRC tissues (Figure 6B). The cBioPortal online database includes 446 KIRC samples with complete genetic alteration data, 89 (20%) of them with alteration in the five genes, and the mRNA upregulation was the most frequent alteration type (Figure 6C). Besides, the KM survival curves of the five SARS-CoV-2–related genes were generated by the Kaplan–Meier plotter (splitting the patients by the best cut-off) to determine their relationship with OS. The results showed that a higher expression of ACADM (p = 7.7e-13) in KIRC predicted better OS, but instead the KIRC patients with higher expression levels of CENPF (p = 4.9e-09), KDELC1 (p = 0.0013), PLOD2 (p = 1.8e-06) and TRMT1 (p = 5.2e-10) exhibited worse OS (Figure 6D). Moreover, it seems that these genes and their somatic copy number alteration were associated with the infiltration of B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells in KIRC (Supplementary Figure S3). The correlation between the expression of the five genes is presented in Supplementary Figure S4; there was the greatest positive correlation between KDELC1 and PLOD2 (r = 0.39) and the greatest negative correlation between ACADM and TRMT1 (r = −0.45).

To further confirm the prognostic potential of the SARS-CoV-2–related gene signature in KIRC with different clinicopathological factors, we conducted series of subgroup survival analyses based on age, gender, grade, and stage in the entire TCGA-KIRC cohort. Survival curves revealed that patients in the high-risk group had a worse outcome compared to those in the low-risk group among different age (≤65 and >65), gender (male and female), grade (G1, G2, G3, and G4), and stage (I, II, III, and IV) subgroups (Figures 7A–H). The clinical characteristics of each patient in the entire TCGA-KIRC cohort and their corresponding risk score are summarized in Figure 7I, including the five model gene expression profiles. We also found that dead patients were more likely to have higher risk scores than alive patients (p < 0.001), suggesting that a higher risk score predicts a poorer prognosis (Figure 7J). For males, the risk scores were significantly higher than females (p < 0.001) (Figure 7K). Furthermore, the risk score was particularly relevant to the tumor malignant phenotype; patients with a higher tumor stage and historical grade usually had higher risk scores, which demonstrated that the prognostic signature could predict the progression of KIRC (Figures 7L,M).

In order to investigate the distribution patterns of patients in high- and low-risk groups, according to the expression of genes in the signature, we performed PCA analyses in the training cohort, testing cohort, and entire KIRC cohort. The results in the three cohorts all showed that patients from different risk groups were distributed in two distinct directions (Figures 8A–C).

Viral infection is associated with inflammation; SARS-CoV-2 infection could induce severe inflammatory responses and cytokine storms, contributing to serious events. To further illuminate the underlying biological processes related to the risk score signature, GSEA was implemented between high- and low-risk groups in the entire TCGA-KIRC cohort. We found some antiviral defense and inflammatory/immune-related processes were active in patients with high risks (Figures 8D–I), such as positive regulation of the defense response to the virus by the host (NES = 2.08, FDR = 0.003), regulation of the acute inflammatory response (NES = 2.17, FDR = 0.002), cytokine production involved in the inflammatory response (NES = 1.62, FDR = 0.03), lymphocyte migration (NES = 2.96, FDR < 0.001), T-cell differentiation involved in the immune response (NES = 2.77, FDR < 0.001), and regulation of T-cell–mediated immunity (NES = 2.02, FDR = 0.004). Taken together, the 5-gene risk score signature based on SARS-CoV-2–related genes could reflect the immune status and deduct the immune landscapes of KIRC patients, and high-risk patients were more likely to have a strong inflammatory response to SARS-CoV-2.

To obtain comprehensive insights into the infiltration of tumor-infiltrating immune cells in KIRC, the CIBERSORT algorithm combined with the LM22 leukocyte signature matrix was utilized to calculate the proportions of 22 immune cell types in high- and low-risk KIRC patients from the entire cohort. In total, 404 patients from the entire TCGA-KIRC cohort met the inclusion criteria (p < 0.05), and their composition of 22 immune cells is summarized in Figure 9A. Besides, the proportions of different tumor-infiltrating immune cell subsets were moderately to weakly correlated (Figure 9B). Between the two groups, the fractions of immune cells in KIRC varied significantly. The KIRC patients in the high-risk group had significantly higher fractions of macrophages M0 (p < 0.001), T follicular helper cells (p < 0.001), regulatory T cells (p < 0.001), dendritic cells activated (p < 0.01), and plasma cells (p < 0.05), and lower fractions of dendritic cells resting (p < 0.001), mast cells resting (p < 0.001), macrophages M2 (p < 0.01), macrophages M1 (p < 0.05), and monocytes (p < 0.05) than patients with low risk (Figure 9C).

In addition, we further assessed the tumor microenvironment of KIRC based on the ESTIMATE method and found that patients in the high-risk group had significantly higher stromal scores (p = 0.0056) and immune scores (p < 0.001) than those in the low-risk group (Figure 10A). At present, the immune checkpoint blockade therapy has become an important antitumor strategy by reversing the immunosuppressive state. Then, the relationships between the expression levels of several key immune checkpoint molecules (CTLA-4, LAG-3, TIGIT, TIM3, and PDCD1) and risk scores were assessed. We observed that the risk scores were significantly positively associated with the expression of CTLA-4, LAG-3, TIGIT, and PDCD1 (p < 0.05) (Figure 10B), and the expression levels of CTLA-4, LAG-3, TIGIT, and PDCD1 were significantly higher in the high-risk KIRC patients than those in the low-risk patients (p < 0.001) (Figure 10C).

In addition, to further extend the clinical applicability and availability of the prognostic model, a predictive nomogram for 1-, 3-, and 5-year OS combined with age, gender, grade, and risk scores was developed in the training set, where predictors were included in the nomogram (Figure 11A). Calibration plots demonstrated that the nomogram had a similar performance than the ideal model for predicting 1-, 3-, and 5-year OS probability, and there was a good agreement between the predicted and actual observed survival (Figure 11B). The KIRC patients in the training set (n = 252) were segregated into high-nomogram score (n = 126) and low-nomogram score (n = 126) groups determined by the median nomogram score, and Kaplan–Meier analysis indicated that the patients with low-nomogram scores tend to have a significantly worse survival (p < 0.001) (Figure 11C). The AUCs of the nomogram in the 1-, 3-, and 5-year ROC curves were 0.880, 0.804, and 0.781, which outperformed the 5-gene signature (Figure 11F; Supplementary Table S5).

Additionally, we further validated the predictive nomogram in the testing set (n = 249) and the entire KIRC cohort (n = 501), and similar results were observed in Kaplan–Meier analysis and ROC curve analysis (Figures 11D,E,G–H; Supplementary Table S5). Altogether, these results indicated that the risk score signature in combination with other clinical parameters might have higher predictive power for OS.

In order to further determine the reliability and robustness of our prognostic signature and nomogram, we used the E-MTAB-1980 dataset for external validation. Table 2 displayed the clinical characteristics of KIRC patients in the E-MTAB-1980 dataset. Risk scores of these patients were calculated using the previous following formula. Similar to our results above, the K-M survival analysis indicated the OS of patients in the high-risk group was significantly shorter than that in the low-risk group (p < 0.001) (Figure 12A). The AUC values of the 5-gene prognostic signature for 1-, 3-, and 5-year OS reached 0.798, 0.783, and 0.807, respectively (Figure 12B). Moreover, the univariate and multivariate Cox regression analyses also demonstrated that the 5-gene signature could predict KIRC patients’ prognosis independently (HR = 2.659, 95% CI = 1.375–5.141, p = 0.004) in Figures 12C,D. We then validated the nomogram in the E-MTAB-1980 dataset and found that KIRC patients with high nomogram scores had poorer prognosis than those with low nomogram scores (p < 0.001) (Figure 12E). The AUC of the nomogram for 5-year OS was 0.875, which exceeded the 5-gene signature and other clinical features (Figure 12F). These results in the E-MTAB-1980 dataset further confirmed the reliability and robustness of the 5-gene prognostic signature and nomogram.