- 1Department of Nephrology, Fujian Provincial Hospital, Shengli Clinical Medical College of Fujian Medical University, Fuzhou, China
- 2Department of Radiation Oncology, Fujian Cancer Hospital, Fujian Medical University Cancer Hospital, Fuzhou, China
Methylation is one of the most extensive modifications of biological macromolecules and affects cell-fate determination, development, aging, and cancer. Several methylation modifications, including 5-methylcytosine and N6-methyladenosine, play an essential role in many cancers. However, little is known about the relationship between methylation and the prognosis of clear cell renal cell carcinoma (ccRCC). Here, we established a methylation-regulating genes prognostic signature (MRGPS) to predict the prognoses of ccRCC patients. We obtained ccRCC samples from The Cancer Genome Atlas and identified methylation-regulatingd genes (MRGs) from the Gene Set Enrichment Analysis database. We also determined differentially expressed genes (DEGs) and performed cluster analysis to identify candidate genes. Subsequently, we established and validated an MRGPS to predict the overall survival of ccRCC patients. This was also verified in 15 ccRCC samples collected from the Fujian Provincial Hospital via quantitative real-time transcription (qRT-PCR). While 95 MRGs were differentially expressed (DEGs1) between tumor and normal tissues, 17 MRGs were differentially expressed (DEGs2) between cluster 1 and 2. Notably, 13 genes common among DEGs1 and DEGs2 were identified as hub genes. In fact, we established three genes (NOP2, NSUN6, and TET2) to be an MRGPS based on their multivariate Cox regression analysis coefficients (p < 0.05). A receiver operating characteristic curve analysis confirmed this MRGPS to have a good prognostic performance. Moreover, the MRGPS was associated with characteristics of the tumor immune microenvironment and responses to inhibitor checkpoint inhibitors. Data from “IMvigor 210” demonstrated that patients with a low MRGPS would benefit more from atelozumab (p < 0.05). Furthermore, a multivariate analysis revealed that MRGPS was an independent risk factor associated with ccRCC prognosis (p < 0.05). Notably, a nomogram constructed by combining with clinical characteristics (age, grade, stage, and MRGPS risk score) to predict the overall survival of a ccRCC patient had a favorable predictive value. Eventually, our qRT-PCR results showed that tumor tissues had higher NOP2 and NSUN6 expression levels and lower TET2 expression than normal tissues of ccRCC samples. While the proposed MRGPS comprising NOP2, NSUN6, and TET2 can be an alternative prognostic biomarker for ccRCC patients, it is a promising index for personalized ICI treatments against ccRCC.
Introduction
Clear cell renal cell carcinoma (ccRCC) is one of the most lethal malignancies of the genitourinary tract, accounting for 70–80% of renal cell carcinoma patients (Zhao et al., 2018). Despite substantial advances in the diagnosis and treatment of ccRCC, long-term prognosis remains far from satisfactory (Siegel et al., 2017). Approximately 20–30% of patients initially present with metastasis (Reiter et al., 2015), indicating that the current screening index for ccRCC is inadequate; thus, it is necessary to immediately identify an aggressive diagnostic marker for ccRCC. In addition, approximately 30–40% of patients with localized ccRCC relapse or exhibit metastasis within 2 years of undergoing radical surgeries (Miao et al., 2018). This implies that the ccRCC patient population is greatly heterogeneous and highlights the inaccuracies in the existing staging system integrated with clinicopathological characteristics.
Interestingly, ccRCC is a highly immunogenic tumor characterized by an abundance of suppressed immune cells (Díaz-Montero et al., 2020). A randomized phase II study has demonstrated that immune checkpoint inhibitor (ICI) monotherapy exhibits non-inferiority efficacy to sunitinib (Mcdermott et al., 2018). However, a CheckMate-214 trial (Cella et al., 2019; Albiges et al., 2020) has revealed that nivolumab combined with iplimumab has positive outcomes compared with sunitinib. Thus, this combination has been approved by the United States Food and Drug Administration as a frontline therapeutic approach for ccRCC patients with intermediate severity. Nonetheless, the objective response rates (ORRs) of avelumab, pembrolizumab, and nivolumab are 16, 36, and 17%, respectively (Tzeng et al., 2021), whereas that of avelumab combined with nivolumab is 42% (Cella et al., 2019). Additionally, continuing treatment with nivolumab has been found to be associated with reduced tumor burden in approximately 50% of patients (Hellmann et al., 2018). Hence, an aggressive biomarker, except PD-1/PD-L1, tumor mutation burden (TMB), and microsatellite status, is urgently warranted in ICI management for ccRCC.
Methylation is one of the most abundant modifications that is widespread across all biological processes. It involves an alkylation reaction, wherein a methyl group replaces a hydrogen atom (Michalak et al., 2019). Methyltransferases, also called “writers,” use the methyl donor S-adenosylmethionine to catalyze methylation; “writers” cooperate with dedicated “erasers” (demethylases) and methyl “readers” (Dawson and Kouzarides, 2012). Genomic studies have demonstrated that hypo- and/or hyper-methylation occur in various enzymes and can result in loss of histone modification (Michalak et al., 2019). Few examples include mutations in metabolic enzymes that regulate histone and DNA demethylation and somatic mutations in core histone genes (You and Jones, 2012). In fact, previous studies have demonstrated that aberrant changes in DNA or RNA methylation can be prospectively utilized in the diagnosis, prognosis, and individualized treatment of various cancers, including ccRCC (Fang et al., 2020; Zhang et al., 2020; Li et al., 2021). Therefore, we systematically analyzed the transcriptomic data of ccRCC patient tissues to identify methylation-regulating genes (MRGs) and accurately predict the prognoses and guide the ICI management of ccRCC patients.
Materials and Methods
Patients and Datasets
We retrieved 359 human MRGs from the Gene Set Enrichment Analysis (GSEA) database (https://www.gsea-msigdb.org/gsea/index.Jsp; Supplementary Table S1) (Subramanian et al., 2005). Moreover, we obtained RNA sequencing (RNA-Seq) expression profile dataset of 537 ccRCC patients and 72 corresponsonding normal samples from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) (Tomczak et al., 2015). The clinicopathological characteristics and survival data of these patients was also retrieved from TCGA. The RNA-seq profiles and clinical data of“IMvigor 210” cohort were obtained from http://research-pub.gene.com/IMvigor210CoreBiologies/.
Furthermore, 15 frozen, surgically resected tumor specimens were acquired from patients pathologically diagnosed with ccRCC at the Fujian Provincial Hospital (FPH) between December 2018 and December 2020. Additionally, we validated the immunohistochemical staining of prognostic genes using The Human Protein Atlas (HPA) database (http://www.proteinatlas.org/) (Uhlén et al., 2015). This study was approved by the ethics committee of the FPH.
Identification of Methylation-Regulating Hub Genes
Based on the RNA-seq data of the ccRCC samples (537 tumors vs 72 normal samples) obtained from TCGA, we analyzed the differentially expressed genes (DEGs1) between tumor and normal tissues. We also functionally explored the biological properties of MRGs in the TCGA ccRCC patients by clustering ccRCC patients into different clusters using the “ConsensusClusterPluspackage” (Wilkerson and Hayes, 2010) (http://www.bioconductor.org/; 1,000 iterations and resampling rate of 80%). The cumulative distribution function (CDF) and delta area were considered to determine the optimal number of groups (k). Subsequently, we identified DEGs between the different clusters (DEGs2) and defined the hub genes as genes common to both DEGs1 and DEGs2.
Construction and Validation of Methylation-Regulating Genes Prognostic Signature
We divided TCGA patients into training and validation cohorts at a ratio of 3:7 (11 samples were deleted because their OS was 0 or unknown), and prognostically significant hub genes (p < 0.05) were screened by univariate Cox regression analysis. In fact, these candidate genes were used to establish a methylation-regulating genes prognostic signature (MRGPS) via multivariate Cox regression analysis. The risk score for each patient was determined using the following formula:
Thereafter, the patients were classified into low-risk and high-risk groups based on the median risk score. We determined the prognostic ability of the MRGPS in the training cohort by generating Kaplan–Meier survival curves and receiving operating characteristic (ROC) curves using the R packages “survminer” and “survivalROC”. The prognostic performance of this MRGPS was further tested in the testing cohort in the same manner as mentioned above.
Functional Analysis
We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to analyze the main function using the “clusterProfiler" (Yu et al., 2012) R package and visualized it using the “Treemap” (Liu et al., 2021) and “ggplot2” packages. Additionally, GSEA was performed to understand the biological processes prevalent in the different subgroups using the “clusterProfiler” R package. Predefined gene sets were identified in the GSEA using the GO Biological Process; 5,000 permutations were performed to determine the p values of these gene sets. Significant pathways were defined as having a p value of <0.05 and a false discovery rate (FDR) of <0.05 (Powers et al., 2018).
Immune Score and Immunotherapy Benefits Analyses
We conducted a single sample GSEA (ssGSEA) analysis, where (Bustin and Mueller, 2005)in we analyzed 20 immune cells of 537 ccRCC samples based on the expression profile of a single sample; we used the “gsva” R package to perform this analysis (Hänzelmann et al., 2013). The ESTIMATE algorithm (i.e., the “estimate” R package) was used to calculate the immune score of each patient. Subsequently, we assessed the immune score difference between the two cluster subgroups. A semi-quantitative analysis of 22 immune cell types in the two MRGPS groups was performed using CIBERSORT via the “cibersort” R package (Chen et al., 2018). Moreover, we calculated tumor immune dysfunction and exclusion (TIDE) and microsatellite instability (MSI) scores from the website of http://tide.dfci.harvar.edu to assess the potential efficacy of ICIs in the two MRGPS subgroups (Fu et al., 2020). We also compared the somatic mutations between the two MRGPS subgroups by obtaining the TMB, i.e., the total number of somatic mutations.
Predicting the Benefits of MRGPS for Tyrosine Kinase Inhibitors
Since VEGFR-targeted therapy remains the first line of treatment for ccRCC, we explored the sensitivity of TKIs, such as sunitinib, sorafenib, pazopanib, and axitinib, stratified by MRGPS. The sensitivity of each TKI was evaluated by IC50 calculation using the “pRRophetic” package (Geeleher et al., 2014), and the corresponding data were obtained from the Genomics of Drug Sensitivity in Cancer database (Yang et al., 2013).
Development of Risk Prediction Model
Furthermore, we conducted a multivariate Cox analysis to evaluate whether the signature-based risk score was independent of other clinical characteristics. The testing cohort was used to further test the performance of the signature in the same manner mentioned above. Thereafter, we generated a nomogram consisting of the current MRGPS and clinical characteristics with p < 0.1. This helped predict the 1-, 3- and 5-years overall survival (OS) of the TCGA ccRCC patients using the “rms” package. Additionally, we evaluated this nomogram using the calibration curve, ROC curve, and decision-making curve (DCA).
Quantitative Reverse Transcription PCR
Relative quantitation of the 15 paired mRNAs was determined by quantitative reverse transcription polymerase chain reaction (qRT-PCR; SuperScript IV Reverse Transcriptase 18090010; Thermo Fisher, United States). The amplification reactions were performed as described previously (Bustin and Mueller, 2005). NSUN6-specific primers were: forward primer, 5′-ATCTGCGTCCGTTTCACC-3′ and reverse primer, 5′-GCTTCCACCACACCTCATC-3'. NOP2-specific primers were: forward primer, 5′-GGGCACAGACACACAAACA-3′ and reverse primer, 5′-GAACGGATGGGAGACACAG-3'. TET2-specific primers were: forward primer, 5′-CACAACCATCCCAGAGTTCA-3′ and reverse primer, 5′-ACTTCCTCCAGTCCCATTTG-3'. Human β-actin-specific primers were: forward primer, 5′-GAAGAGCTACGAGCTGCCTGA-3′ and reverse primer 5′- CAGACAGCACTGTGTTGGCG-3'. Data analysis was performed using the ΔΔCT method.
Statistical Analyses
Distributed data were compared by performing the Student’s t-test and Wilcoxon test, whereas proportion differences were calculated by the chi-square test. Additionally, component analysis in subgroups were compared by the Fisher’s test. While survival differences between different groups were assessed via the log-rank test, prognostic factors were identified by the Cox regression analyses. All statistical analyses were performed using RStudio version 4.0.3, and two-sided p < 0.05 was considered as statistically significant.
Results
Identification of Methylation-Regulating Hub Genes
The entire analytical process of this study is presented in Figure 1, and the clinical characteristics of the ccRCC patients in the TCGA and FPH cohorts are listed in Table 1. Additionally, Figure 2A presents all 95 DEGs among tumor and normal tissues (DEGs1: FDR<0.05 and |log2FC|>0.5), including 51 upregulated and 44 downregulated genes; the top 50 DEGs are presented in a heatmap (Figure 2B). The GO analysis of DEGs1 revealed that methylation-relate biological process (BP), cellular component (CC), and molecular function (MF) were enriched in tumor tissues (Figure 2C).
FIGURE 2. Identification of methylation-regulating Hub Genes. (A) Volcano plot demonstrates DEGs1. (B) Heatmap demonstrates the top 50 DEGs1. (C) The Gene Ontology (GO) analysis of DEGs1 in ccRCC. (D) Consensus clustering cumulative distribution function (CDF) for k = 2 to 9. (E) The consensus heatmap showed the ccRCC patients was divided into two distinct clusters when k = 2. (F) Overall Survival (OS) analysis of different clusters in the TCGA dataset. (G) Volcano plot demonstrates the DEGs2. (H) Heatmap demonstrates the DEGs2. (I) Venn diagram demonstrates the intersect between DEGs1 and DEGs2.
Subsequently, we performed consensus clustering to explore the molecular characteristics between different MRG expression samples. We observed a relative change in the CDF of the consensus cluster from k = 2 to k = 9 (Figure 2D); the delta area under the CDF curve from k = 2 to 9 is depicted in Supplementary Figure S1H. The corresponding heatmap presents the results of this consensus from k = 2 to 9 (k = 2, Figure 2E; k = 3–9, Supplementary Figure S1A–G). The criteria for deciding the cluster number was determined by a relatively high consistency and a low variation coefficient and an appreciable increase in the area under the CDF curve. Thus, after comprehensive consideration, we chose k = 2 as the optimal cut-off for the clusters number.
A significant difference in the OS was observed between patients of clusters 1 and 2 (p < 0.001, Figure 2F). To determine which genes contributed to this difference in prognosis, we first identified 17 genes as DEGs between cluster 1 and cluster 2 (DEGs2: FDR<0.05, |log2FC|>0.5; Figure 2G); these genes are also represented in a heatmap (Figure 2H). Eventually, 13 overlapping genes between DEGs1 and DEGs2 were identified as the hub genes (Figure 2I, Supplementary Table S2).
We further evaluated the molecular characteristics of the different clusters by conducting immune-related analyses between clusters 1 and 2. The ssGSEA demonstrated that cluster 2 had a high abundance of approximately all immune cell types compared to cluster 1 (p < 0.05, Figure 3A). Moreover, the tumor microenvironment estimate scores, including the stromal, immune, and total scores, were higher in cluster 2 than those in cluster 1 (p < 0.05, Figure 3B). Notably, immune-related signaling pathways were enriched in cluster 2, as determined by the GSEA (Figure 3C and Supplementary Table S3). Furthermore, we obtained the enrichment score of each negative immune-related signaling pathway in each ccRCC sample by performing a gene set variation analysis (GSVA). Consequently, we observed significant survival differences between high- and low-GSVA scores regarding the negative regulation of adaptive immune response, negative regulation of immune response, and negative regulation of leukocyte-mediated immunity (p < 0.05, Figures 3D–F); however, this did not hold true for the negative regulation of natural killer cell-mediated immunity (p > 0.05, Figure 3G).
FIGURE 3. Immune cell infiltration analysis and GSEA analysis between different clusters. (A) Estimated abundance of 20 immune cells using ssGSEA. (B) Tumor microenvironment (TME) estimate score in different clusters. (C) GSEA delineation of the biological pathways which enrich in cluster 2 using the gene set “c5. go.bp.v7.4. symbols”. Overall Survival (OS) analysis in different GSVA score of (D) negative regulation of adaptive immune response, (E) negative regulation of immune response, and (F) negative regulation of leukocyte-mediated immunity, (G) negative regulation of natural killer cell-mediated immunity in TCGA-ccRCC patients. Significant statistical differences between the two clusters were assessed using the Wilcoxon test (ns, p > 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001).
Construction and Validation of the Methylation-Regulating Genes Prognostic Signature
In the training cohort, we screened prognosis-associated seven hub genes by univariate Cox regression analysis (Figure 4A). Then, a multivariate Cox regression analysis was conducted to screen the optimal model and was depicted in Figure 4B; based on their regression coefficients, three MRGS (NOP2, NSUN6, and TET2) were identified to form an MRGPS. The MRGPS score of each patient was calculated according to the following formula: Risk score = [NOP2 expression*(0.656940513)] + [NSUN6 expression*(0.911107243)] + [TET2 expression*(-1.180533124)]. Considering the median score as the cut-off value, patients in the training cohort were divided into low- and high-risk groups; these patients had apparent survival differences (p < 0.001, Figure 4C). The corresponding risk scores and survival statuses are presented in Figure 4D. The ROC curves demonstrated the excellent predictive capability of the current MRGPS with 1-, 3-, and 5-years AUCs being 0.798, 0.750, and 0.768, respectively (Figure 4E). Likewise, the advantages of the current MRGPS were observed in the validation (Figures 4F–H) and the whole cohorts (Figures 4I,J).
FIGURE 4. Construction and validation of the MRGPS. (A) Forrest plot of the univariate Cox regression analysis in the training cohort. (B) Forrest plot of the multivariate Cox regression analysis in the training cohort. Kaplan-Meier analysis, risk score analysis and ROC curve of the MRGPS inthe training cohort (C–E), validation cohort (F–H), and whole TCGA cohort (I–K).
In addition, relationships between clinicopathological characteristics and risk scores were further explored. As shown in Supplementary Figure S2A, differences were observed regarding the age (age ≤65 years, age >65 years), differentiation (G1, G2, G3, G4), T stage (T1, T2, T3, T4), M stage (M0, M1), and cancer stage (I, II, III, IV). Furthermore, Kaplan–Meier survival curves showed that the high-risk patients had worse prognoses than the low-risk patients in the following attributes: age ≤65 years, age >65 years, male sex, female sex, G1-2, G3-4, T1-2, T3-4, M0, M1, stage I–II, and stage III–IV (p < 0.05, Supplementary Figure S2B–M).
Immune Analyses and Immunotherapy
We further explored the immune microenvironment characteristics of patients belonging to the different risk subgroups by conducting immune cell infiltration and immune function analysis on the TCGA cohort patients. We observed significantly decreased number of naive B cells, memory B cells, plasma cells, CD4+T cells, CD4+T memory cells, gamma T cells, resting NK cells, M0/M1/M2 and resting dendritic cells, activated dendritic cells, and resting mast cells and also observed increased number of CD8+T cells and regulatory T cells in the high-risk group (p < 0.05, Figure 5A). Relative expression levels of MHC molecules and co-stimulatory molecules and adhesion factors, such as CD40, CD58, HLA-A, HLA-B, HLA-C, HLA-DMA, HLA-DOB, HLA-DPB1, and HLA-F, were all higher in the high-risk group than those in the low-risk group (p < 0.05, Figure 5B). Importantly, the expression levels of immune checkpoint proteins, such as PDCD1, CTLA4, TBX2, TNF, LAG3, CD8A, IFNG, and GZMB were all significantly higher in the high-risk group than those in the low-risk group (p < 0.05, Figure 5C).
FIGURE 5. Immune-related analysis between different MRGPS subgroup. (A) The proportions of TME cells in different MRGPS subgroups. (B) Relative expression of MHC molecules, co-stimulatory molecules, and adhesion factors. (C) Association of MRGPS with immune checkpoint molecules. (D) TIDE, (E) MSI, (F) TMB score in different MRGPS subgroups. (G) Distribution of immune response to ICIs therapy in different MRGPS subgroups in IMvigor patients. Significant statistical differences between the two subgroups were assessed using the Wilcoxon test (ns, p > 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001).
Furthermore, patients in the high-risk group were found to have a higher TIDE score, lower MSI score, and higher TMB than those in the low-risk group (p < 0.05, Figures 5D–F). This suggested that low-risk patients may benefit more from immunotherapy compared to high-risk ones according to the current MRGPS. We further validated this observation using the “IMvigor 210” dataset containing clinical information and RNA-seq data of metastatic urothelial cancer patients who were treated with the ICI atezolizumab (PD-L1 inhibitor). Remarkably, data from the 298 IMvigor patients also validated the clinical utility of the current MRGPS in response to atezolizumab (p < 0.05, Figure 5G).
Potential Biological Pathway Analysis of Methylation-Regulating Genes Prognostic Signature
We further determined the potential biological pathways prevalent in different risk group by performing a KEGG analysis on the DEGs among high- and low-risk groups (FDR<0.05 and |log2FC|>0.5). We observed that the “PI3K−Akt signaling pathway”, “mTOR signaling pathway”, “Ras signaling pathway” and other carcinogenesis-related pathways were enriched in the high-risk group (Figure 6A). Furthermore, the GSEA revealed that the high-risk group had higher enrichment score for the PI3K-AKT and mTOR pathways compared to the low-risk group (Figures 6B,C). These results revealed that the MRGPS possibly promotes cancer development by activating these pathways.
FIGURE 6. KEGG and GSEA analysis of MRGPS. (A) KEGG analysis of the DEGs between high- and low-risk groups. (B) PI3K−Akt signaling pathway and (C) mTOR signaling pathway were identified in the high-risk group.
VEGF Family Expressions and TKI Sensitivity
As the VEGF family was an important molecular target, we compared their expression levels in high- and low-risk groups. Consequently, no significant differences in VEGFA expression were observed between the two groups (p > 0.05, Figure 7A). However, VEGFB and VEGFD expression levels were significantly upregulated in the high-risk group (p < 0.05, Figures 7B,D), whereas that of VEGFC was significantly downregulated in the high-risk group (p < 0.05, Figure 7C). Further analysis revealed that sunitinib had lower IC50 the higher-risk group than that in the low-risk groups (p < 0.05, Figure 7E). In contrast, pazopanib (p < 0.05, Figure 7G), but not sorafenib and axitinib (both p > 0.05, Figures 7F,H), had lower IC50 in the low-risk group than that in the high-risk group. These results indicated that different risk groups had varying susceptibilities for different targeted drugs.
FIGURE 7. VEGF family expressions and the sensibility of TKI inhibitors in different MRGPS subgroups. (A–D) VEGF family expressions in different MRGPS subgroups. Drug susceptibility analysis between different MRGPS subgroups about (E) Sunitinib, (F) Sorafenib, (G) Pazopanib, (H) Axitinib. Wilcoxon test (ns, p > 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001).
Methylation-Regulating Genes Prognostic Signature-Based Nomogram Construction
Furthermore, we discovered that the current risk score was an independent risk factor for OS in the training, validation, and whole cohorts (p < 0.05, Table 2). Subsequently, we developed a nomogram based on age, differentiation grade, stage, and MRGPS risk score to further predict the OS of ccRCC patients belonging to the TCGA cohort (p < 0.1, Figure 8A). We observed good calibrations regarding the predicted vs observed 1-, 3-, and 5-years OS of the patients (Figure 8B). Moreover, the ROC curves exhibited better predictive capability in the current nomogram to predict the 1-, 3-, and 5-years OS than the MRGPS and risk scores published by Wang et al. (2021), Chen et al. (2021), and Zheng et al. (2021) (Figures 8C–E). Additionally, DCA analysis revealed the superiority of the current nomogram over MRGPS and the published risk scores in predicting the 1-, 3-, and 5-years OS (Figures 8F–H).
FIGURE 8. Construction and verification of nomogram. (A) The prognostic nomogram constructed based on the risk score of MRGPS and clinicopathological parameters predicted the survival rate of TCGA-ccRCC patients at 1-, 3-, and 5-years. (B) Calibration curves showed the concordance between predicted and observed 1-, 3-, and 5-years survival rates. AUCs of the nomogram, MRGPS and other signatures in ROC analysis were calculated at (C) 1-, (D) 3-, and (E) 5-years OS time in TCGA-ccRCC cohort. Decision curve analyses (DCA) for nomogram, MRGPS and other signatures at (F) 1-, (G) 3-, and (H) 5-years to assess clinical utility in TCGA-ccRCC cohort.
Validation Using Quantitative Real-Time Transcription-PCR and Human Protein Atlas Datasets
Our qRT-PCR analysis revealed elevated expression levels of NOP2 and NSUN6, but decreased expression of TET2 were in the tumor tissues compared to those in the paired normal tissues of 15 ccRCC samples obtained from FPH (p < 0.05, Figures 9A–C). The results of HPA database demonstrated that the expression levels of both NOP2 and NSUN6 were higher in the ccRCC tissues than those in the normal tissues; however, the expression of TET2 was significantly lower in the ccRCC tissues than that in the normal tissue (Figures 9D–F).
FIGURE 9. Validation using qPCR and HPA datasets. (A) NOP2, (B) NSUN6, and (C) TET2 mRNA expression measured by qRT-PCR. Validation of the differences in expression of (D) NOP2, (E) NSUN6, and (F) TET2 between renal cancer and normal renal tissue at the translational level with data from the HPA database.
Discussion
Global and local changes in DNA/RNA/histone methylation are seminal features of malignant tumor cells (Michalak et al., 2019). In the current study, we identified three MRGs (NOP2, NSUN6, and TET2) from TCGA data and established an MRGPS for the prognoses of ccRCC patients. This MRGPS exhibited excellent calibration and discrimination. In addition, we validated the three candidate genes in 15 paired ccRCC samples obtained from FPH by qRT-PCR. Furthermore, the current risk score was correlated with tumor immune microenvironment characteristics and could be used as a potential biomarker of ccRCC response to ICIs.
Of note, ccRCC is a highly heterogeneous malignancy (Jonasch et al., 2021). The existing prognosis models that incorporate clinicopathological characteristics, such as the AJCC staging system and the Mayo Clinic stage and necrosis score, have improved prognosis capacity (Parker et al., 2017). However, owing to the complex molecular mechanism of ccRCC, clinical parameters alone are inadequate for predicting the prognoses of ccRCC patients. Interestingly, chromatin methylations, such as m5C and m6A, play a fundamental role in the ccRCC carcinogenesis (Angulo et al., 2021). Nonetheless, comprehensive exploration of chromatin methylation in ccRCC is still lacking. In this study, we established a novel MRGPS using data from TCGA ccRCC patients; this MRGPs improved the prognoses of ccRCC patients with a C-index as high as 0.798 at 1-year OS. In addition, close links were identified between the clinical and pathological characteristics of ccRCC and MRGPS: age, sex, differentiation, and tumor node metastasis (TNM) stage. Furthermore, a MRGPS-incorporating nomogram demonstrated a higher prognostic capacity and clinical utility than published risk scores.
Among the three MRGs identified, NOP2 and NSUN6 were prognostic risk factors, whereas TET2 was a prognostic protective factor. Notably, NOP2 and NSUN6 are key members of the NOP2/Sun domain family and possess S-adenosyl-L-methionine-dependent methyltransferase activity (Frye and Blanco, 2016). NOP2 is upregulated in various cancers, including lung adenocarcinoma, breast cancer, and prostate cancer, and it is associated with tumor aggressiveness (Ma et al., 2017). Deficiency of NSUN6-mediated methylation can downregulate transcription and translation. While NSUN6 expression is highest in the testis and lowest in the blood, it is heterogeneous in different tumors. However, it is downregulated in tumors originating from tissues that have high NSUN6 expression, such as the testis, thyroid, and ovaries (Selmi et al., 2021). In contrast, it is upregulated in tumors originating from tissues that have low NSUN6 expression, such as that in hematologic tumor and kidney cancer. Moreover, NSUN6 is associated with prognosis of various cancers, including pancreatic cancer (Yang et al., 2021) and hepatocellular carcinoma (Wang et al., 2018). It also plays an important role in bone metastasis (Li et al., 2017). On the other hand, TET2 mutations have been widely identified various myeloid malignancies. In fact, TET2 inactivation leads to polyhematopoietic abnormalities in mice, which is a recurrent event in human lymphoma formation (Ferrone et al., 2020). Notably, TET2 dysfunction mutations are generally associated with DNA hypermethylation, tumor progression, and poor patient outcomes (Cimmino et al., 2017). However, NOP2, NSUN6, and TET2 have been rarely studied in ccRCC. In the present study, qRT-PCR data from 15 paired FPH ccRCC samples revealed that while NOP2 and NSUN6 were upregulated, TET2 was downregulated in tumor tissues compared with those in normal tissues. In summary, NOP2, NSUN6, and TET2 were identified as prognostic biomarkers for ccRCC; however, additional in vitro and in vivo research is needed to validate these findings.
The potential mechanisms of MRGPS regulating ccRCC prognosis deserved further study. In the present study, we found via KEGG analysis and GSEA that the PI3K-AKT and mTOR signaling pathways were highly enriched in the high-risk subgroup. The PI3K signaling pathway facilitates several essential cellular functions, such as cell proliferation, growth, migration, metabolism, and survival (Fruman and Rommel, 2014). In a large cohort of 419 primary ccRCC patients, aberrantly expressed components of the PI3K signaling cascade (e.g., PTEN, PI3K, p-AKT, mTOR, p-mTOR, p-S6, and p-4EBP1 proteins) exhibited aggressive pathological features and caused adverse survival (Darwish et al., 2013). Therefore, we hypothesized that poor prognoses of patients with a high MRGPS might be because of activation of the PI3K-AKT and mTOR pathways; nonetheless, this hypothesis requires further exploration.
Since ccRCC is a highly vascular tumor, the levels of angiogenic factors, including VEGF, are correlated with its prognosis (Choueiri and Kaelin, 2020). Inhibition of VEGFR generally causes vascular normalization, thereby activating anti-tumor immunity (Hsieh et al., 2017). Until 2017, the multikinase inhibitors sunitinib and pazopanib that primarily target VEGFR formed the frontline treatment for ccRCC (Powles et al., 2021). The median progress free survival (PFS), OS and ORR for sunitinib and pazopanib are 8.4 and 9.5 months, 28.4 months and 29.3 months, and 25 and 31%, respectively (George et al., 2021). The complexity of the VEGFR family may possibly be responsible for the inconsistent results. In this study, we found that the current MRGPS could be used as an alternative to VEGFRs. Remarkably, it had a positively correlation with VEGFB/D and negative correlation with VEGFC; however, it did not have a correlation with VEGFA. Furthermore, we also found that while sunitinib had a lower IC50, pazopanib had a higher IC50 in high-risk patients than those in low-risk patients, according to the current MRGPS. This highlighted the response divergence between sunitinib and pazopanib and clarified personalized TKI treatment for ccRCC patients.
Although ccRCC patients have a typically suppressed immune status, they are highly abundant in immune cells (Şenbabaoğlu et al., 2016). In this study, we revealed that the current MRGPS was correlated with tumor-infiltrating lymphocytes: high MRGPS was associated with increased number of CD8+T cell number, activated NK cell, follicular helper cells T cells, and regulatory T cells. In contrast, low MRGPS was associated with decreased number of naïve B cells, resting memory CD4+ T cells, monocytes, macrophages M2, and resting mast cells, as previously reported (Díaz-Montero et al., 2020). In addition, the current MRGPS was also associated with co-stimulatory molecules, such as CD40 and CD58. Immune checkpoints are cell surface receptors expressed on immune cells, and their inhibition causes immune activation. In the present study, we found that high MRGPS was associated with PDCD1 and CTLA4 expression levels (both p < 0.001). Further analysis revealed that low MRGPS was correlated with lower TIDE score and higher MSI score than high MRGPS (both p < 0.05), indicating that patients with low MRGPS would benefit more from ICIs. Importantly, this finding was validated in 298 IMvigor patients receiving atezolizumab. Therefore, MRGPS is a promising biomarker for predicting the response of ccRCC patients towards ICIs.
Nonetheless, our study has several limitations. First, both the data from TCGA and FPH are retrospective; therefore, the risk score needs to be verified in prospective cohorts. Second, merely incorporating MRGPS to build a prognostic model is inadequate, regardless of its importance. Third, samples from FPH were too few, and the results need to be validated in more samples. Finally, the associations of the current MRGPS with tumor mutations, tumor immune microenvironment, and TKI and ICI responses require further validation in vitro and in vivo.
In conclusion, the current MRGPS consisting of NOP2, NSUN6, and TET2 is a potential alternative prognostic biomarker for ccRCC patients and is also be a promising index for personalized ICI treatments in ccRCC.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: Expression profile data are available in The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Immunohistochemical staining of genes are available in The Human Protein Atlas (HPA) (http://www.proteinatlas.org/). And the experimental data of qRT-RCR can be downloaded from the supplementary file.
Author Contributions
LZ, ZS, LW and FH were involved in the study concept and design and data analysis and interpretation and in drafting the manuscript. All authors contributed to the article and approved the submitted version.
Funding
Supported by Natural Science Foundation of Fujian Province of China (No. 2019J01185) and key Clinical Specialty Discipline Construction Program of Fujian, P.R.C (No. 2017739).
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.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We sincerely thank the researchers and study participants for their contributions towards this study.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.832803/full#supplementary-material
Supplementary Figure S1 | The ccRCC patients from TCGA was divided into different clusters when k=3-9 (A–G). Relative change in area under CDF curve for k=2 to 9 (H).
Supplementary Figure S2 | The correlation analysis with clinical characteristic. (A) Bargraph demonstrates the proportion of patients in the high- and low- risk groups with different clinicopathological characteristics. (B-M) Kaplan-Meier analysis of subgroup patients based on risk scores.
References
Albiges, L., Tannir, N. M., Burotto, M., Mcdermott, D., Plimack, E. R., Barthélémy, P., et al. (2020). Nivolumab Plus Ipilimumab versus Sunitinib for First-Line Treatment of Advanced Renal Cell Carcinoma: Extended 4-year Follow-Up of the Phase III CheckMate 214 Trial. ESMO Open 5, e001079. doi:10.1136/esmoopen-2020-001079
Angulo, J. C., Manini, C., López, J. I., Pueyo, A., Colás, B., and Ropero, S. (2021). The Role of Epigenetics in the Progression of Clear Cell Renal Cell Carcinoma and the Basis for Future Epigenetic Treatments. Cancers (Basel) 13, 2071. doi:10.3390/cancers13092071
Bustin, S. A., and Mueller, R. (2005). Real-time Reverse Transcription PCR (qRT-PCR) and its Potential Use in Clinical Diagnosis. Clin. Sci. (Lond) 109, 365–379. doi:10.1042/cs20050086
Cella, D., Grünwald, V., Escudier, B., Hammers, H. J., George, S., Nathan, P., et al. (2019). Patient-reported Outcomes of Patients with Advanced Renal Cell Carcinoma Treated with Nivolumab Plus Ipilimumab versus Sunitinib (CheckMate 214): a Randomised, Phase 3 Trial. Lancet Oncol. 20, 297–310. doi:10.1016/s1470-2045(18)30778-2
Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12
Chen, H., Pan, Y., Jin, X., and Chen, G. (2021). Identification of a Four Hypoxia-Associated Long Non-coding RNA Signature and Establishment of a Nomogram Predicting Prognosis of Clear Cell Renal Cell Carcinoma. Front. Oncol. 11, 346. doi:10.3389/fonc.2021.713346
Choueiri, T. K., and Kaelin, W. G. (2020). Targeting the HIF2-VEGF axis in Renal Cell Carcinoma. Nat. Med. 26, 1519–1530. doi:10.1038/s41591-020-1093-z
Cimmino, L., Dolgalev, I., Wang, Y., Yoshimi, A., Martin, G. H., Wang, J., et al. (2017). Restoration of TET2 Function Blocks Aberrant Self-Renewal and Leukemia Progression. Cell 170, 1079–1095. e1020. doi:10.1016/j.cell.2017.07.032
Darwish, O. M., Kapur, P., Youssef, R. F., Bagrodia, A., Belsante, M., Alhalabi, F., et al. (2013). Cumulative Number of Altered Biomarkers in Mammalian Target of Rapamycin Pathway Is an Independent Predictor of Outcome in Patients with clear Cell Renal Cell Carcinoma. Urology 81, 581–586. doi:10.1016/j.urology.2012.11.030
Dawson, M. A., and Kouzarides, T. (2012). Cancer Epigenetics: from Mechanism to Therapy. Cell 150, 12–27. doi:10.1016/j.cell.2012.06.013
Díaz-Montero, C. M., Rini, B. I., and Finke, J. H. (2020). The Immunology of Renal Cell Carcinoma. Nat. Rev. Nephrol. 16, 721–735. doi:10.1038/s41581-020-0316-3
Fang, J., Hu, M., Sun, Y., Zhou, S., and Li, H. (2020). Expression Profile Analysis of m6A RNA Methylation Regulators Indicates They Are Immune Signature Associated and Can Predict Survival in Kidney Renal Cell Carcinoma. DNA Cel Biol 39, 2194–2211. doi:10.1089/dna.2020.5767
Ferrone, C. K., Blydt-Hansen, M., and Rauh, M. J. (2020). Age-Associated TET2 Mutations: Common Drivers of Myeloid Dysfunction, Cancer and Cardiovascular Disease. Int. J. Mol. Sci. 21, 626. doi:10.3390/ijms21020626
Fruman, D. A., and Rommel, C. (2014). PI3K and Cancer: Lessons, Challenges and Opportunities. Nat. Rev. Drug Discov. 13, 140–156. doi:10.1038/nrd4204
Frye, M., and Blanco, S. (2016). Post-transcriptional Modifications in Development and Stem Cells. Development 143, 3871–3881. doi:10.1242/dev.136556
Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., et al. (2020). Large-scale Public Data Reuse to Model Immunotherapy Response and Resistance. Genome Med. 12, 21. doi:10.1186/s13073-020-0721-z
Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS One 9, e107468. doi:10.1371/journal.pone.0107468
George, D. J., Lee, C. H., and Heng, D. (2021). New Approaches to First-Line Treatment of Advanced Renal Cell Carcinoma. Ther. Adv. Med. Oncol. 13, 17588359211034708. doi:10.1177/17588359211034708
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Hellmann, M. D., Callahan, M. K., Awad, M. M., Calvo, E., Ascierto, P. A., Atmaca, A., et al. (2018). Tumor Mutational Burden and Efficacy of Nivolumab Monotherapy and in Combination with Ipilimumab in Small-Cell Lung Cancer. Cancer Cell 33, 853–861. e854. doi:10.1016/j.ccell.2018.04.001
Hsieh, J. J., Purdue, M. P., Signoretti, S., Swanton, C., Albiges, L., Schmidinger, M., et al. (2017). Renal Cell Carcinoma. Nat. Rev. Dis. Primers 3, 17009. doi:10.1038/nrdp.2017.9
Jonasch, E., Walker, C. L., and Rathmell, W. K. (2021). Clear Cell Renal Cell Carcinoma Ontogeny and Mechanisms of Lethality. Nat. Rev. Nephrol. 17, 245–261. doi:10.1038/s41581-020-00359-2
Li, C., Wang, S., Xing, Z., Lin, A., Liang, K., Song, J., et al. (2017). A ROR1-HER3-lncRNA Signalling axis Modulates the Hippo-YAP Pathway to Regulate Bone Metastasis. Nat. Cel Biol 19, 106–119. doi:10.1038/ncb3464
Li, H., Jiang, H., Huang, Z., Chen, Z., and Chen, N. (2021). Prognostic Value of an m5C RNA Methylation Regulator-Related Signature for Clear Cell Renal Cell Carcinoma. Cancer Manag. Res. 13, 6673–6687. doi:10.2147/cmar.s323072
Liu, L., Chandrashekar, P., Zeng, B., Sanderford, M. D., Kumar, S., and Gibson, G. (2021). TreeMap: a Structured Approach to fine Mapping of eQTL Variants. Bioinformatics 37, 1125–1134. doi:10.1093/bioinformatics/btaa927
Ma, P., Pan, Y., Li, W., Sun, C., Liu, J., Xu, T., et al. (2017). Extracellular Vesicles-Mediated Noncoding RNAs Transfer in Cancer. J. Hematol. Oncol. 10, 57. doi:10.1186/s13045-017-0426-y
Mcdermott, D. F., Huseni, M. A., Atkins, M. B., Motzer, R. J., Rini, B. I., Escudier, B., et al. (2018). Clinical Activity and Molecular Correlates of Response to Atezolizumab Alone or in Combination with Bevacizumab versus Sunitinib in Renal Cell Carcinoma. Nat. Med. 24, 749–757. doi:10.1038/s41591-018-0053-3
Miao, D., Margolis, C. A., Gao, W., Voss, M. H., Li, W., Martini, D. J., et al. (2018). Genomic Correlates of Response to Immune Checkpoint Therapies in clear Cell Renal Cell Carcinoma. Science 359, 801–806. doi:10.1126/science.aan5951
Michalak, E. M., Burr, M. L., Bannister, A. J., and Dawson, M. A. (2019). The Roles of DNA, RNA and Histone Methylation in Ageing and Cancer. Nat. Rev. Mol. Cel Biol 20, 573–589. doi:10.1038/s41580-019-0143-1
Parker, W. P., Cheville, J. C., Frank, I., Zaid, H. B., Lohse, C. M., Boorjian, S. A., et al. (2017). Application of the Stage, Size, Grade, and Necrosis (SSIGN) Score for Clear Cell Renal Cell Carcinoma in Contemporary Patients. Eur. Urol. 71, 665–673. doi:10.1016/j.eururo.2016.05.034
Powers, R. K., Goodspeed, A., Pielke-Lombardo, H., Tan, A.-C., and Costello, J. C. (2018). GSEA-InContext: Identifying Novel and Common Patterns in Expression Experiments. Bioinformatics 34, i555–i564. doi:10.1093/bioinformatics/bty271
Powles, T., Albiges, L., Bex, A., Grünwald, V., Porta, C., Procopio, G., et al. (2021). ESMO Clinical Practice Guideline Update on the Use of Immunotherapy in Early Stage and Advanced Renal Cell Carcinoma. Ann. Oncol.. 32:1511–1519. doi:10.1016/j.annonc.2021.09.014
Reiter, M. A., Kurosch, M., and Haferkamp, A. (2015). Renal Cell Carcinoma: Drug Therapy and Prognostic Models. Urologe A 54, 735–738. 735-746; quiz. doi:10.1007/s00120-015-3845-9
Selmi, T., Hussain, S., Dietmann, S., Heiß, M., Borland, K., Flad, S., et al. (2021). Sequence- and Structure-specific Cytosine-5 mRNA Methylation by NSUN6. Nucleic Acids Res. 49, 1006–1022. doi:10.1093/nar/gkaa1193
Şenbabaoğlu, Y., Gejman, R. S., Winer, A. G., Liu, M., Van Allen, E. M., De Velasco, G., et al. (2016). Tumor Immune Microenvironment Characterization in clear Cell Renal Cell Carcinoma Identifies Prognostic and Immunotherapeutically Relevant Messenger RNA Signatures. Genome Biol. 17, 231. doi:10.1186/s13059-016-1092-z
Siegel, R. L., Miller, K. D., and Jemal, A. (2017). Cancer Statistics, 2017. CA: A Cancer J. Clinicians 67, 7–30. doi:10.3322/caac.21387
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102, 15545–15550. doi:10.1073/pnas.0506580102
Tomczak, K., Czerwińska, P., and Wiznerowicz, M. (2015). The Cancer Genome Atlas (TCGA): an Immeasurable Source of Knowledge. Contemp. Oncol. (Pozn) 19, A68–A77. doi:10.5114/wo.2014.47136
Tzeng, A., Tzeng, T. H., and Ornstein, M. C. (2021). Treatment-free Survival after Discontinuation of Immune Checkpoint Inhibitors in Metastatic Renal Cell Carcinoma: a Systematic Review and Meta-Analysis. J. Immunother. Cancer 9, e003473. doi:10.1136/jitc-2021-003473
Uhlén, M., Fagerberg, L., Hallström, B. M., Lindskog, C., Oksvold, P., Mardinoglu, A., et al. (2015). Tissue-based Map of the Human Proteome. Science 347, 1260419. doi:10.1126/science.1260419
Wang, Z.-L., Li, B., Luo, Y.-X., Lin, Q., Liu, S.-R., Zhang, X.-Q., et al. (2018). Comprehensive Genomic Characterization of RNA-Binding Proteins across Human Cancers. Cel Rep. 22, 286–298. doi:10.1016/j.celrep.2017.12.035
Wang, B., Zhao, H., Ni, S., and Ding, B. (2021). Integrated Analysis of the Roles of RNA Binding Proteins and Their Prognostic Value in Clear Cell Renal Cell Carcinoma. J. Healthc. Eng. 2021, 5568411. doi:10.1155/2021/5568411
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics 26, 1572–1573. doi:10.1093/bioinformatics/btq170
Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of Drug Sensitivity in Cancer (GDSC): a Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111
Yang, R., Liang, X., Wang, H., Guo, M., Shen, H., Shi, Y., et al. (2021). The RNA Methyltransferase NSUN6 Suppresses Pancreatic Cancer Development by Regulating Cell Proliferation. EBioMedicine 63, 103195. doi:10.1016/j.ebiom.2020.103195
You, J. S., and Jones, P. A. (2012). Cancer Genetics and Epigenetics: Two Sides of the Same coin? Cancer Cell 22, 9–20. doi:10.1016/j.ccr.2012.06.008
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Zhang, D., Wang, Y., and Hu, X. (2020). Identification and Comprehensive Validation of a DNA Methylation-Driven Gene-Based Prognostic Model for Clear Cell Renal Cell Carcinoma. DNA Cel Biol. 39, 1799–1812. doi:10.1089/dna.2020.5601
Zhao, H., Cao, Y., Wang, Y., Zhang, L., Chen, C., Wang, Y., et al. (2018). Dynamic Prognostic Model for Kidney Renal clear Cell Carcinoma (KIRC) Patients by Combining Clinical and Genetic Information. Sci. Rep. 8, 17613. doi:10.1038/s41598-018-35981-5
Keywords: methylation, clear cell renal cell carcinoma, quantitative real-time transcription, immune checkpoint inhibitor, prognosis, risk signature
Citation: Zhang L, Su Z, Hong F and Wang L (2022) Identification of a Methylation-Regulating Genes Prognostic Signature to Predict the Prognosis and Aid Immunotherapy of Clear Cell Renal Cell Carcinoma. Front. Cell Dev. Biol. 10:832803. doi: 10.3389/fcell.2022.832803
Received: 10 December 2021; Accepted: 04 February 2022;
Published: 02 March 2022.
Edited by:
Biaoru Li, Augusta University, United StatesReviewed by:
Shukui Wang, Nanjing Medical University, ChinaBorbala Mifsud, Hamad bin Khalifa University, Qatar
Copyright © 2022 Zhang, Su, Hong and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fuyuan Hong, hongdoc@163.com; Lei Wang, wangleiy001@126.com