key: cord-0739725-lcku0ive
authors: Huang, Yongbiao; Chen, Sheng; Xiao, Lingyan; Qin, Wan; Li, Long; Wang, Yali; Ma, Li; Yuan, Xianglin
title: A Novel Prognostic Signature for Survival Prediction and Immune Implication Based on SARS-CoV-2–Related Genes in Kidney Renal Clear Cell Carcinoma
date: 2022-01-24
journal: Front Bioeng Biotechnol
DOI: 10.3389/fbioe.2021.744659
sha: 3063129216bbdbbca33e25c87a7e4d78b2af65f6
doc_id: 739725
cord_uid: lcku0ive

Kidney renal clear cell carcinoma (KIRC) is a common aggressive malignancy of the urinary system. COVID-19, a highly infectious and severe disease caused by SARS-CoV-2, has become a significant challenge for global public health. Cancer patients have been reported to be more vulnerable to SARS-CoV-2 infection and have a higher risk for serious complications than the general population. However, the correlation between KIRC and COVID-19 remains incompletely elucidated. In this study, we comprehensively investigated the expression and prognostic significance of 333 SARS-CoV-2 infection–related genes in KIRC using the TCGA dataset and identified 31 SARS-CoV-2–related differently expressed genes between KIRC and normal renal tissues. Based on these genes, we constructed and validated a 5-gene prognostic signature (including ACADM, CENPF, KDELC1, PLOD2, and TRMT1) to distinguish low- and high-risk KIRC patients of poor survival in TCGA and E-MTAB-1980 cohorts. Gene set enrichment analysis (GSEA) showed that some inflammatory/immune-related pathways were significantly enriched in the high-risk group. The ESTIMATE analysis indicated that patients in the high-risk group had higher stromal and immune cell scores, therefore lower tumor purity. Moreover, they presented higher proportions of macrophages M0, regulatory T cells (Tregs), and T follicular helper cells and higher expression of immune checkpoints CTLA-4, LAG-3, TIGIT, and PDCD1 than low-risk patients. Besides, we also developed a nomogram to expand clinical applicability, which exhibits excellent predictive accuracy for survival. In conclusion, we identified a novel prognostic signature and nomogram based on SARS-CoV-2–related genes as reliable prognostic predictors for KIRC patients and provided potential therapeutic targets for KIRC and COVID-19.

Renal cell carcinoma (RCC) has been considered as the third most common malignancy of the urinary system in adults, with high morbidity and mortality (Cohen and McGovern, 2005; Rini et al., 2009 ). More than 300,000 new cases are diagnosed, and over 100,000 patients die from it per year worldwide; it has gradually become a global public health burden (Siegel et al., 2020) . Kidney renal clear cell carcinoma (KIRC) is the most frequent histopathologic subtype of RCC, accounting for approximately 75% of all primary RCCs (Jonasch et al., 2014; Moch et al., 2016) . Despite recent improvements in multiple therapeutic methods, including radical resection, chemoradiotherapy, targeted therapy, and immunotherapy, the prognosis for most KIRC patients remains poor because of the high risk of local recurrence and distant metastasis (Barata and Rini, 2017; Choueiri and Motzer, 2017) . The recurrence rate of KIRC after initial curative nephrectomy varies from 20 to 30% within 5 years, and the 5year overall survival rate of metastatic KIRC patients is less than 10%. Up to 30% of KIRC patients already had metastases at the first time of diagnosis and missed the optimal treatment opportunity (Cohen and McGovern, 2005; Rini et al., 2009) . Thus, it is imperative to explore novel reliable biomarkers for early diagnosis and prognosis prediction and provide prognostic indicators and treatment targets for KIRC.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a new coronavirus infecting humans, was first identified in Wuhan, China, in December 2019 and spread rapidly around the world and caused a global coronavirus disease 2019 pandemic within several months (Tian et al., 2020; . As of December 26, 2021, more than 279 million confirmed cases of SARS-CoV-2 infection and 5.39 million COVID-19-related deaths have been reported worldwide, based on the latest statistics from Johns Hopkins University (https://coronavirus.jhu.edu/map.html).

Epidemiological studies so far have demonstrated that patients diagnosed with cancer have a higher risk of SARS-CoV-2 infection and resulted in serious complications and poor prognosis, especially in patients undergoing anticancer treatments (chemotherapy, surgery, or radiotherapy) and patients with hematologic malignancies (Dai et al., 2020; Liang et al., 2020; Wang and Zhang, 2020; Yu et al., 2020) .

A recent study revealed that the SARS-CoV2 receptors or auxiliary proteins ACE2, ANPEP, ENPEP, and DPP4 are highly expressed in renal tumor, especially in KIRC. Furthermore, these SARS-CoV2 receptors are closely correlated with immune infiltrates and immune response in KIRC (Mihalopoulos et al., 2020; Tripathi et al., 2020) . However, besides the several aforementioned SARS-CoV-2-related proteins, more than 300 SARS-CoV-2-interacting human proteins have been identified (Gordon et al., 2020) , and the correlations between their expression, prognosis, and immune cell infiltration in KIRC tissues have not been elucidated. Thus, our study was carried out to explore the expression levels, prognostic value, and immune perspective of these SARS-CoV-2-related proteins in KIRC and further construct a SARS-CoV-2-related prognostic signature for KIRC. This study revealed the important roles of SARS-CoV-2-interacting proteins in KIRC and illustrated the susceptibility of KIRC patients to SARS-CoV-2 infection and the underlying mechanism of COVID-19 in KIRC patients.

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 n i (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 timedependent 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 .

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 flowchart of our study was presented in Figure 1 .

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 deadenylationdependent 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 ).

Validating the 5-Gene Risk Score Signature in the Testing Set and Entire TCGA-KIRC Cohort

To further examine the predictive ability of this risk score model, we also calculated the risk scores of the KIRC patients in the Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 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 cutoff) 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 Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 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 highand low-risk groups in the entire TCGA-KIRC cohort. We found some antiviral defense and inflammatory/immunerelated 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 tumorinfiltrating 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 Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 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 Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 8 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 lownomogram 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 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.

COVID-19 pandemic is becoming a serious global public health problem; the older patients with basic diseases such as hypertension, diabetes, chronic kidney disease, and cancer were considered at an increased risk of COVID-19 infection (Ma et al., 2020; Salzberger et al., 2021) . This susceptibility may Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 also be related to the aberrant expression of some key SARS-CoV-2-interacting proteins, especially ACE2 Ashraf et al., 2021; Rezaei et al., 2021) . SARS-CoV-2 utilizes ACE2 as an entry receptor to infect cells, and an increased ACE2 expression may result in an elevated SARS-CoV-2 infection risk (Cao et al., 2020; Wu et al., 2021) . Some studies have found that ACE2 is aberrantly expressed in some cancers, such as lung adenocarcinoma, renal carcinoma, and stomach adenocarcinoma, and it is correlated with prognosis and immune infiltrates (Hoang et al., 2020; Kong et al., 2020; Tripathi et al., 2020; Chen L. et al., 2021; He et al., 2021) . Moreover, the other key SARS-CoV-2 infection-related host cell receptors or entry related-proteins, TMPRSS2, CTSL/B, USP13, HSPA5, Furin, and ADAM17 have recently been reported to play important roles in SARS-CoV-2 infection Ibrahim et al., 2020; Zipeto et al., 2020; Tang et al., 2021) . They were also observed to be aberrantly expressed in several cancers; dysregulation of their expression in cancer patients' tissues should affect the susceptibility and severity of SARS-CoV-2 infection (Li H. et al., 2020; Ravaioli et al., 2020; Süt, 2020; Fu et al., 2021; Zhou and Gao, 2021) . These findings suggest that SARS-CoV-2 infection-related genes may play essential roles in virus infection, and they also participate in development of many cancers. To date, over 300 SARS-CoV-2 infection-related genes Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 have been identified, but their roles in KIRC remains unclear. It is important to evaluate their expression in KIRC tissues, for predicting the KIRC patients' susceptibility to SARS-CoV-2 infection and outcome. Therefore, the goal of our study was to develop a robust SARS-CoV-2-related gene signature, which can predict the clinical outcome for KIRC patients.

In the present study, we systematically analyzed the expression of 333 SARS-CoV-2-related genes in the TCGA-KIRC cohort and found 31 SARS-CoV-2-related DEGs between KIRC and normal tissues. Moreover, we further investigated the association between these SARS-CoV-2-related DEGs and overall survival of KIRC patients by using univariate Cox regression, multiple stepwise Cox regression, and K-M survival analyses and generated a novel prognostic model comprising five SARS-CoV-2-related genes. In the training set, the median risk score was used as the cut-off to categorize KIRC patients into high-risk or low-risk groups, and the high-risk patients had unfavorable prognosis. Multivariate analyses demonstrated that the risk score signature was an independent prognostic factor for OS prediction of KIRC patients. The ROC curves exhibited that the AUC values of the SARS-CoV-2-related gene signature were greater than 0.7 at 1, 3, and 5 years, indicating that the 5gene signature had a good predictive performance to screen out the KIRC patients with poor prognosis. Additionally, we further validated the risk score model in the testing set, entire TCGA-KIRC cohort, and E-MTAB-1980 dataset and observed the similar results and equally good performance. Stratification analysis suggested that the risk score signature had wide applicability for predicting prognosis of KIRC patients with distinct molecular features and clinicopathological characteristics. The 5-gene signature was able to accurately distinguish the KIRC patients who had a high risk with poor overall survival in different clinical subgroups. All these results indicated that our risk score model was a robust indicator for predicting KIRC patients' prognosis.

The prognostic signature which we constructed consisted of five SARS-CoV-2-related genes (ACADM, CENPF, KDELC1, PLOD2, and TRMT1). Of the five genes, CENPF, KDELC1, PLOD2, and TRMT1 were upregulated in KIRC tissues and were correlated with poor survival and higher susceptibility to SARS-CoV-2. In contrast, ACADM was downregulated and act as Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 a protective factor for KIRC patients; the higher ACADM expression may contribute to better survival and lower likelihood of SARS-CoV-2 infection. Centromere protein F (CENPF), a gene located in chromosome 1q41, encodes a microtubule-binding protein that is associated with the cell cycle and regulate chromosome segregation; it begins to accumulate during the S-phase and reaches maximal expression levels in the G2/M phase cells (Shahid et al., 2018; Sun et al., 2019) . CENPF has been reported to be upregulated in various types of malignant tumors and serve as a prognostic indicator for multiple cancers, such as prostate cancer, breast cancer, pancreatic cancer, lung adenocarcinoma, and renal cancer; the high expression of CENPF is associated with worse prognosis and metastasis in these cancers . Some studies have found that CENPF and FOXM1 are major regulators in prostate cancer development; they synergistically promote prostate cancer malignant progression and metastasis (Aytes et al., 2014) . Moreover, CENPF also promotes breast cancer progression and bone metastasis by activating the AKT/mTOR signaling pathway . Knockdown of CENPF could inhibit the progression of lung adenocarcinoma by suppressing the ERβ2/5 pathway (Hexiao et al., 2021) . KDELC1, also called POGLUT2, is a protein O-glucosyltransferase which participates in O-glucose modification and modulates Notch trafficking and signaling (Takeuchi et al., 2018; Mehboob and Lang, 2021) . Currently, few research studies are available regarding the role of KDELC1 in cancer. In this research, we first uncovered that KDELC1 was upregulated in KIRC and associated with poor survival. PLOD2 is one of the members of the PLOD family, encodes for lysyl hydroxylases 2which mediates stabilized collagen cross-link formation in the extracellular matrix. Many lines of evidence suggested that PLOD2 is a tumor metastasis-promoting gene; it promotes tumor aggressive metastasis and invasion through mediating cross-links of collagen (Yamauchi and Sricholpech, 2012; Chen et al., 2015) . In recent years, increasing research has revealed the pivotal role of PLOD2 in various cancer types, particularly in breast cancer, hepatocellular carcinoma, bladder cancer, sarcomas, and renal cell carcinoma (Du et al., 2017; Qi and Xu, 2018) . In these cancer tissues, PLOD2 is significantly overexpressed, and the high PLOD2 expression was significantly associated with tumor metastasis and shorter survival time. Obviously, PLOD2 is an independent factor of poor outcomes and could serve as a prognostic biomarker for patients with these cancer types (Du et al., 2017) . Several previous studies have shown that the expression PLOD2 could be regulated by HIF-1α (Gilkes et al., 2013) . Besides, in a tumor microenvironment, PLOD2 could modulate the collagen-cross-linking activity of cancerassociated fibroblasts and promote tumor cell migration and invasion (Chen et al., 2015; Qi and Xu, 2018) . TRMT1 encodes tRNA methyltransferase 1, an RNA methyltransferase responsible for N 2 , N 2 -dimethylguanosine (m2,2G) formation in cytosolic and mitochondrial tRNAs (Dewe et al., 2017; Jonkhout et al., 2021) . TRMT1 has been considered as the cause of autosomal-recessive intellectual disability , but its role in cancers has not yet been reported. Herein, our study has reported the prognostic value of TRMT1 in KIRC for the first time. ACADM encodes medium-chain acyl-CoA dehydrogenase which is Frontiers in Bioengineering and Biotechnology | www.frontiersin.org January 2022 | Volume 9 | Article 744659 15 involved in fatty acid oxidation and lipid metabolism (Smith et al., 2010) . Inhibition of ACADM could alter hepatocellular carcinoma cell lipid metabolism and promote cell lipid accumulation, ultimately drive hepatocarcinogenesis . Additionally, a recent study showed that ACADM was highly expressed in ferret models after SARS-CoV-2 infection and played a crucial role in SARS-CoV-2 infection progression . In summary, two of the prognostic model genes (CENPF and PLOD2) have been reported to be associated with KIRC patients' prognosis in previous studies, while the remaining three genes (KDELC1, TRMT1, and ACADM) were first reported to be correlated with the prognosis of KIRC in the current study.

Evidence has suggested that the immune dysregulation may be highly linked to the pathological process of COVID-19. Given the significance of the immune system in antiviral and antitumor responses, we further investigated the immune cell infiltration status between KIRC patients with low and high risks based on the prognostic signature. The ESTIMATE analysis showed that patients in the high-risk groups were identified to have significantly higher stromal scores and immune scores and lower tumor purities. Recent studies have revealed that major stromal components within the tumor microenvironment can not only favor tumor growth and metastasis but it could also affect the antitumor immune response, leading to unfavorable prognosis (Quail and Joyce, 2013; Turley et al., 2015) . Our results were in agreement with these previous findings that high stromal infiltration and low tumor purity levels correlated with worse prognosis in many cancers (Zhang et al., 2017) . Furthermore, the CIBERSORT algorithm was used to evaluate the abundance of 22 immune cell types in KIRC, and the higher proportions of macrophages M0, regulatory T cells (Tregs), and T follicular helper cells were observed in high-risk patients than low-risk patients. Tumor-associated macrophages and Treg cells are immunosuppressive cells that inhibit the antitumor immune response and secrete various immunosuppressive cytokines, promoting the tumor immune escape (Dunn et al., 2002; Sica et al., 2006) . Moreover, we noted that the high-risk KIRC patients had significantly higher levels of the CTLA-4, LAG-3, TIGIT, and PDCD1 expressions than the low-risk patients, and high expression levels of immune checkpoints are effective predictors for responses to immune checkpoint inhibitors (Pardoll, 2012; Qin et al., 2019) . Notably, Treg cells also expressed immune checkpoint molecules, including CTLA-4 and PDCD1 (Romano et al., 2015) . The above results indicated that the worse survival of the KIRC patients with high-risk scores is likely due to the higher immune checkpoint expression levels and more potent immunosuppressive tumor microenvironment. Besides, these results also suggested that patients in the high-risk group will benefit more from the immune checkpoint blockade therapy than those in the low-risk group, thereby contributing to a better prognosis.

Additionally, we further established a novel predictive nomogram by integrating our risk signature with several clinical characteristics with superior prediction performance. The AUCs of this nomogram for 1-, 3-, and 5-year survival are all greater than 0.750 in the training set, testing set, entire TCGA-KIRC cohort, and E-MTAB-1980 dataset, which indicated our nomogram and prognostic signature have excellent predictive ability for both short-and long-term follow-up patients. We also found the AUC values of this nomogram exceeded that of the 5-gene signature both in predicting 1-, 3-, and 5-year OS, indicating combination of the 5-gene signature with other prognosis-related clinical parameters would further improve the predictive power for survival and expand the clinical practicability.

To our knowledge, our study for the first time proposed a novel prognostic signature and nomogram for KIRC based upon SARS-CoV-2-related genes and provided some new insights into the relationship between KIRC and SARS-CoV-2 infection. However, some limitations of our study still existed. First, the study was based on retrospective cohorts; our results needed to be further confirmed by future prospective studies. Second, the underlying biological functions and molecular mechanisms of the five model genes in KIRC and SARS-CoV-2 infection have not been examined in this study and remain to be further studied in the future.

In summary, we developed a novel prognostic signature and nomogram based on five SARS-CoV-2-related genes for the first time, which exhibited good performance both in the TCGA cohort and external validation dataset. This prognostic model was confirmed to have independent prognostic significance for KIRC patients and may provide some new potential therapeutic targets for KIRC patients and protects them from SARS-CoV-2 infection. 

SARS-CoV-2, ACE2 Expression, and Systemic Organ Invasion

Cross-species Regulatory Network Analysis Identifies a Synergistic Interaction between FOXM1 and CENPF that Drives Prostate Cancer Malignancy

Nomograms in Oncology: More Than Meets the Eye

Treatment of Renal Cell Carcinoma: Current Status and Future Directions

Comparative Genetic Analysis of the Novel Coronavirus

Receptor ACE2 in Different Populations

Centromere Protein F Is Identified as a Novel Therapeutic Target by Genomics Profile and Contributing to the Progression of Pancreatic Cancer

Single Cell Transcriptome Revealed SARS-CoV-2 Entry Genes Enriched in colon Tissues and Associated with Coronavirus Infection and Cytokine Production

Lung Adenocarcinoma Patients Have Higher Risk of SARS-CoV-2 Infection

Lysyl Hydroxylase 2 Induces a Collagen Cross-Link Switch in Tumor Stroma

Systemic Therapy for Metastatic Renal-Cell Carcinoma

Renal-cell Carcinoma

Patients with Cancer Appear More Vulnerable to SARS-COV-2: a Multi-center Study during the COVID-19 Outbreak

TRMT1-Catalyzed tRNA Modifications Are Required for Redox Homeostasis to Ensure Proper Cellular Proliferation and Oxidative Stress Survival

PLOD2 in Cancer Research

Cancer Immunoediting: From Immunosurveillance to Tumor Escape

GeneMANIA Update

Evaluation and Characterization of HSPA5 (GRP78) Expression Profiles in normal Individuals and Cancer Patients with COVID-19

Integrative Analysis of Complex Cancer Genomics and Clinical Profiles Using the cBioPortal

Procollagen Lysyl Hydroxylase 2 Is Essential for Hypoxia-Induced Breast Cancer Metastasis

A SARS-CoV-2-Protein Interaction Map Reveals Drug Targets and Potential Drug Repurposing

Integrated Bioinformatic Analysis of SARS-CoV-2 Infection Related Genes ACE2, BSG and TMPRSS2 in Aerodigestive Cancers

Knockdown of CENPF Inhibits the Progression of Lung Adenocarcinoma Mediated by ERβ2/5 Pathway

Genetic Susceptibility of ACE2 and TMPRSS2 in Six Common Cancers and Possible Impacts on COVID-19

A Novel RNA Binding Protein-Related Prognostic Signature for Hepatocellular Carcinoma

COVID-19 Spike-Host Cell Receptor GRP78 Binding Site Prediction

Renal Cell Carcinoma

Subcellular Relocalization and Nuclear Redistribution of the RNA Methyltransferases TRMT1 and TRMT1L upon Neuronal Activation

Analysis of the Susceptibility of Lung Cancer Patients to SARS-CoV-2 Infection

Genomic, Epigenomic, and Immune Subtype Analysis of CTSL/B and SARS-CoV

Frontiers in Bioengineering and Biotechnology | www.frontiersin.org

TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells

Yin Yang 1 Facilitates Hepatocellular Carcinoma Cell Lipid Metabolism and Tumor Progression by Inhibiting PGC-1β-Induced Fatty Acid Oxidation

Physiological and Pathological Regulation of ACE2, the SARS-CoV-2 Receptor

Cancer Patients in SARS-CoV-2 Infection: A Nationwide Analysis in China

Gene Signatures of SARS-CoV/SARS-CoV-2-Infected Ferret Lungs in Shortand Long-Term Models

Development and Validation of a TP53-Associated Immune Prognostic Model for Hepatocellular Carcinoma

Clinical Characteristics and Prognosis in Cancer Patients with COVID-19: A Single center's Retrospective Study

Structure, Function, and Pathology of Protein O-Glucosyltransferases

Molecular Determinants and Clinical Implications in Renal Cancer

The 2016 WHO Classification of Tumours of the Urinary System and Male Genital Organs-Part a: Renal, Penile, and Testicular Tumours

Pancancer Survival Analysis of Cancer Hallmark

Robust Enumeration of Cell Subsets from Tissue Expression Profiles

The Blockade of Immune Checkpoints in Cancer Immunotherapy

Roles of PLODs in Collagen Synthesis and Cancer Progression

Novel Immune Checkpoint Targets: Moving beyond PD-1 and CTLA-4

Microenvironmental Regulation of Tumor Progression and Metastasis

ACE2 and TMPRSS2 Potential Involvement in Genetic Susceptibility to SARS-COV-2 in Cancer Patients

ACE2: Its Potential Role and Regulation in Severe Acute Respiratory Syndrome and COVID-19

Renal Cell Carcinoma

Ipilimumab-dependent Cell-Mediated Cytotoxicity of Regulatory T Cells Ex Vivo by Nonclassical Monocytes in Melanoma Patients

Epidemiology of SARS-CoV-2

Integrated Molecular Analysis of clear-cell Renal Cell Carcinoma

Modulates Cancer Metabolism by Regulating Pyruvate Kinase M2 Phosphorylation Signaling

Tumour-associated Macrophages Are a Distinct M2 Polarised Population Promoting Tumour Progression: Potential Targets of Anti-cancer Therapy

Cancer Statistics, 2020. CA A

Allelic Diversity in MCAD Deficiency: The Biochemical Classification of 54 Variants Identified during 5years of ACADM Sequencing

Identification of an Immune-Related Long Non-coding RNA Signature and Nomogram as Prognostic Target for Muscle-Invasive Bladder Cancer

Overexpression of CENPF Correlates with Poor Prognosis and Tumor Bone Metastasis in Breast Cancer

Molecular Profiling of Immune Cell-Enriched Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Interacting Protein USP13

Two Novel Protein O-Glucosyltransferases that Modify Sites Distinct from POGLUT1 and Affect Notch Trafficking and Signaling

Proteolytic Activation of SARS-CoV-2 Spike at the S1/S2 Boundary: Potential Role of Proteases beyond Furin

GEPIA2: An Enhanced Web Server for Large-Scale Expression Profiling and Interactive Analysis

Clinical Characteristics and Risk Factors Associated with COVID-19 Disease Severity in Patients with Cancer in Wuhan, China: A Multicentre, Retrospective, Cohort Study

Renal Carcinoma Is Associated with Increased Risk of

Immunological Hallmarks of Stromal Cells in the Tumour Microenvironment

A Novel Coronavirus Outbreak of Global Health Concern

Risk of COVID-19 for Patients with Cancer

Advances in Research on ACE2 as a Receptor for 2019-nCoV

Lysine post-translational Modifications of Collagen

Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data

Frontiers in Bioengineering and Biotechnology | www.frontiersin.org

SARS-CoV-2 Transmission in Patients with Cancer at a Tertiary Care Hospital in Wuhan, china

Tumor Purity as an Underlying Key Factor in Glioma

An Intellectual Disability-associated Missense Variant in TRMT1 Impairs tRNA Modification and Reconstitution of Enzymatic Activity

Pan-Cancer Analysis of FURIN as a

ACE2/ADAM17/TMPRSS2 Interplay May Be the Main Risk Factor for COVID-19

Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

YH and SC designed the study. YH and LX performed data analysis. YW, WQ, LM, and LL revised this manuscript. YH and XY drafted and modified the manuscript. 

We acknowledge the following databases: the HPA database (http://www.proteinatlas.org/), TCGA database (https://portal. gdc.cancer.gov/repository), ArrayExpress database (https:// www.ebi.ac.uk/arrayexpress/), GEPIA database (http://gepia2. cancer-pku.cn/), Kaplan-Meier Plotter database (http://kmplot. com), cBioportal database (http://www.cbioportal.org/), TIMER database (https://cistrome.shinyapps.io/timer/), and Venny2.1 (http://bioinformatics.psb.ugent.be/webtools/ Venn/).

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2021.744659/ full#supplementary-material