A Novel Ferroptosis-Related Pathway for Regulating Immune Checkpoints in Clear Cell Renal Cell Carcinoma.

Ferroptosis is a novel form of cell death and plays a role in various diseases, especially tumors. It has been reported that ferroptosis is involved in the growth and progression of clear cell renal cell carcinoma (ccRCC); however, the specific molecular mechanisms are still unclear. In this study, we constructed a four-gene signature (FeSig) of ferroptosis-related genes via Cox regression analysis. ROC and survival analyses indicated that FeSig had good diagnostic and prognostic value. Further analysis revealed that ferroptosis was associated with tumor immunity in ccRCC. Next, weighted gene co-expression network analysis was performed to identify the potential regulatory mechanisms. Combined with correlation and survival analyses, the TAZ/WNT10B axis was identified as a tumor immune-related regulatory pathway. In conclusion, these findings suggest that ferroptosis is correlated with tumor immunity. The TAZ/WNT10B axis may be a novel biomarker and therapeutic target for immunotherapy in ccRCC.

Introduction

Renal cell carcinoma (RCC) is a common neoplasm that originates from renal tubular epithelial cells and accounts for 2–3% of all malignant tumors in adults (1). In the USA, it is estimated that there will be approximately 70,000 new cases and 14,000 deaths due to this cancer in 2020. The 5-year relative survival is 75.2% (2). Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype and accounts for over 70% of RCCs (1). Due to high resistance to conventional chemoradiotherapy, only limited therapies are currently available (3, 4). At present, immunotherapy, as a novel treatment, is gradually being applied in ccRCC therapy (5). Immune checkpoint inhibitors (ICIs), such as PD-1 and CTLA-4 inhibitors (6, 7), have improved clinical responses and quality of life for some patients. However, some patients are insensitive to ICIs (8). Therefore, it is necessary to further study the regulatory mechanisms of immune checkpoints in ccRCC.

Ferroptosis is a novel form of cell death caused by iron-dependent oxidative damage (9). Due to the failure of glutathione peroxidase (GPX4), a large number of reactive oxygen radicals (ROS) accumulate on membrane lipids (10). In recent years, many studies have verified that ferroptosis plays a vital role in degenerative diseases (i.e., Alzheimer’s and Parkinson’s diseases) (11, 12), ischemia–reperfusion injury (13), and tumors (14). Zhang et al. showed that BAP1 inhibits tumor progression by promoting cellular ferroptosis (15). Furthermore, it has been reported that ferroptosis is closely associated with the therapeutic effect of immunotherapy. Wang et al. verified that CD8+ T cells enhanced the effect of immunotherapy by promoting ferroptosis of tumor cells (16). Interestingly, Lang et al. indicated that immunotherapy sensitizes tumors to radiotherapy by inhibiting SLC7A11 and reducing cystine uptake (17). However, the interaction between ferroptosis and immunotherapy in ccRCC is still unclear. In this study, we constructed a ferroptosis-related gene signature (FeSig) and analyzed the correlation between FeSig and immune checkpoints to identify novel therapeutic targets in ccRCC.

Methods and Materials

Data Download and Study Design

The gene expression data and clinical data were obtained from The Cancer Genome Atlas (TCGA-KIRC, https://www.cancer.gov/tcga) and International Cancer Genome Consortium (ICGC-RECA-EU, https://dcc.icgc.org/) databases. An immunotherapy (immune checkpoint inhibitor, ICI) dataset of clear cell renal cell carcinoma (ccRCC) with clinical data and gene expression data was obtained from the study of Miao et al. (18). In addition, the study process is shown in a flow chart ( ).

Figure 1

Flow chart of study design. The data collection and study process were shown in a flow chart. TCGA, The Cancer Genome Atlas; KIRC, kidney clear cell carcinoma; FRGs, ferroptosis-related genes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KEGG: Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; ESTIMATE, Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data; WGCNA, weighted gene co-expression network analysis; OS, overall survival; DFS, disease-free survival.

Identification of Differentially Expressed Genes and Ferroptosis-Related Genes

The “DESeq2” (19) and “edgeR“ (20) packages were used to screen DEGs with p-value <0.05 and |log2FC| > 1.5. The Wilcoxon test was performed to identify DE-FRGs with a p-value <0.05 and |log2FC| >1.

Protein–Protein Interaction Network Analysis

The DE-FRGs were inputted into the STRING online tool (https://string-db.org) to construct the PPI network. Then, Cytoscape software (21) was used to analyze the PPI network.

Cox Regression Analysis and Construction of a Proportional Hazards Model

The “survival” package (https://CRAN.R-project.org/package=survival) was utilized to perform univariate and multivariate Cox regression analyses. The “glmnet” (22) and “survival” packages were used to perform least absolute shrinkage and selection operator (LASSO) regression analysis. Through integrated Cox analysis, four key FRGs were screened to construct the risk model (FeSig). The risk score (RS) formula was as follows:

Where Coef (i) represents the coefficient, and X(i) represents the expression of selected genes.

Principal Component Analysis and T-Distributed Stochastic Neighbor Embedding Analysis

The “stats” and “limma” packages (23) were used to perform PCA. The “Rtsne” package (https://CRAN.R-project.org/package=Rtsne) was utilized to perform t-SNE analysis. PCA and t-SNE analysis were both used to explore the distribution of different groups.

Survival Analysis and Receiver Operator Characteristic Curve Analysis

According to the median gene expression/risk score, ccRCC patients were divided into a high group and a low group. Then, survival curves of overall survival (OS) and disease-free survival (DFS) were drawn by the “survival” package in R and GraphPad software (version 7.0). A p-value <0.05 was considered statistically significant. Moreover, the “survivalROC” (https://CRAN.R-project.org/package=survivalROC) package was used to generate a time-dependent ROC curve to evaluate the predictive value of the risk model.

Gene Set Enrichment Analysis and Gene Set Variation Analysis

All patients were divided into two groups (high group and low group) according to the median gene expression/risk score. GSEA (24) was performed to discover potential mechanisms and downstream signaling pathways. Moreover, the “GSVA” package (25) was utilized to find differential signaling pathways between the high group and the low group.

Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis

The “clusterProfiler” package (26) was used to conduct KEGG enrichment analysis of DEGs. The results were visualized by the “ggplot2” package (https://CRAN.R-project.org/package=ggplot2) in the R programme. A p-value <0.05 was selected as the cut-off point.

Weighted Gene Co-Expression Network Analysis

The “WGCNA” package (27) was applied to construct the weighted gene co-expression network of DEGs. First, the TCGA-KIRC cohort was evaluated via sample clustering to detect outliers with a height cut-off point = 60. Then, Pearson’s correlation was calculated between each of the gene pairs. Second, a matrix of adjacencies was constructed according to Pearson’s correlation. Then, the adjacencies achieved scale-free topology based on the soft threshold power β ( ). Third, the adjacencies were transformed into a topological overlap matrix (TOM). Genes with a high absolute correlation were clustered into the same module. Finally, combined with clinical traits, we calculated the correlation between the module eigengenes (MEs) and clinical traits to screen the clinically significant modules. After that, the correlation between genes and clinical traits (cor.geneTraitSignificance) and the correlation between genes and MEs (cor.geneModuleMembership) were conducted to identify hub genes. In this study, cor.geneTraitSignificance >0.2 and cor.geneModuleMembership >0.6 were selected as the cut-off criteria.

Figure 7

Identification of key genes of regulating immune checkpoints. (A) Analysis of scale-free fit parameter and mean connectivity for various soft-thresholding power. (B) Dendrogram of DEGs clustered based on a dissimilarity measure (1-TOM). The DEGs were divided into four modules (brown, turquoise, blue, and gray). (C) The correlation between clinical traits and modules. MEturquoise was selected as the hub module. (D–H) Four genes (CPNE7, WNT10B, ADAMTS14, and RUFY4) were regarded as downstream hub genes. (I) The correlation among the FeSig (BID, TAZ, MT1G, and SLC7A11), downstream hub genes (CPNE7, WNT10B, ADAMTS14, and RUFY4), and immune checkpoints (PD-1, CTLA4, LAG3, and TIGIT). BID, TAZ, CPNE7, WNT10B, and RUFY4 were selected for further analysis. DEGs, differentially expressed genes; TOM, topological overlap matrix.

Human Clinical Specimens

A total of 19 pairs of ccRCC samples and adjacent normal samples (4 cm away from the margin of the tumor tissues) were obtained from Wuhan Union Hospital between 2018 and 2020. The study was approved by the Human Research Ethics Committee of Huazhong University of Science and Technology (HUST), and all patients signed the informed consent.

RNA Extraction and qRT-PCR

Total RNA was extracted from ccRCC samples and adjacent normal samples using Trizol Reagent (Sigma, USA). Then, cDNA was synthesized using qPCR RT Kit (Vazyme, China). After that, quantitative real time PCR (qRT-PCR) was performed to amplify cDNA.

Primer sequences were listed as follows:

GAPDH Forward: 5′-GCACCGTCAAGGCTGAGAAC-3′;

GAPDH Reverse: 5′-TGGTGAAGACGCCAGTGGA-3′;

TAZ Forward: 5′-CACCGTGTCCAATCACCAGTC-3′;

TAZ Reverse: 5′-TCCAACGCATCAACTTCAGGT-3′;

WNT10B Forward: 5′-CATCCAGGCACGAATGCGA-3′;

WNT10B Reverse: 5′-CGGTTGTGGGTATCAATGAAGA-3′.

Statistical Analysis

In this study, all data are presented as the mean ± SD. SPSS (version 22.0) and GraphPad Prism (version 7.0) were used to analyze the data. Student’s t-test, Mann–Whitney test, and Pearson’s χ2 test were used to conduct statistical analyses. A p-value <0.05 was considered statistically significant.

Results

Identification of DE-FRGs and Construction of a Proportional Hazards Model

The gene expression data were obtained from the TCGA-KIRC cohort. According to the cut-off criteria, 46 DE-FRGs were identified ( ). The protein–protein interactions of DE-FRGs were shown in . Further analysis indicated that there were 20 DE-FRGs with good prognostic value ( ). The correlation was close and high between these DE-FRGs ( ). Then, LASSO and multivariate Cox regression analyses were applied to construct a prognostic model. A four-gene signature (FeSig) was identified ( ). The risk score = 0.49*expression of BID + 0.70*expression of TAZ + 0.09*expression of MT1G + 0.45*expression of SLC7A11.

Figure 2

Identification of key DE-FRGs via Cox regression analysis. (A) The volcano plot of FRGs. (B) The protein–protein interaction between FRGs. (C) Univariate cox regression analysis revealed that a total of twenty DE-FRGs had prognostic value. (D) The correlation network of DE-FRGs. (E, F) Through LASSO Cox regression, nine genes were selected as key DE-FRGs for further analysis. (G) Four genes (BID, TAZ, MT1G, and SCL7A11) were identified as independent prognostic factors via multivariate Cox regression analysis. DE-FRGs, differentially expressed ferroptosis-related genes; LASSO, least absolute shrinkage and selection operator. ***P < 0.001; *P < 0.05.

The Diagnostic and Prognostic Value of the Four-Gene Signature

PCA and t-SNE analysis proved that FeSig could significantly divide patients into different risk groups ( ). In the TCGA-KIRC cohort, patients were divided into two groups (high-risk group and low-risk group) according to the median risk score ( ). As shown in , patients with a high-risk score had a higher probability of death earlier than those with a low risk score. Further Cox regression analysis revealed that the risk score was an independent prognostic factor ( ). Time-dependent ROC curves indicated that the risk score had good predictive performance in both the TCGA-KIRC ( ) cohort and the ICGC-RECA-EU ( ) cohort. Moreover, survival analysis also uncovered that the high-risk group predicted poor OS ( ). These findings suggested that FeSig had good diagnostic and prognostic value.

Figure 3

The diagnostic and prognostic value of the four-gene signature (FeSig). (A, B) PCA and t-SNE analysis both verified that FeSig could significantly divide patients into different risk groups. (C) The distribution and median value of the risk scores in the TCGA cohort. (D) Patients with high-risk score had poor OS. (E, F) Univariate and multivariate Cox regression analyses indicated that risk score was an independent prognostic factor (HR = 2.48, p < 0.001). (G, H) The risk score had good diagnostic value. (I, J) High-risk group predicted poor OS both in TCGA and ICGC cohort. FeSig, ferroptosis-related gene signature; PCA, principal component analysis; t-SNE, t-distributed stochastic neighbor embedding; TCGA, The Cancer Genome Atlas; OS, overall survival; ICGC, International Cancer Genome Consortium.

FeSig Is Closely Correlated With Immune-Related Pathways

GSEA, GSVA, and KEGG analyses were performed to identify potential downstream signaling pathways of FeSig. As shown in , FeSig was associated with immune cell-related pathways. Analogously, GSVA also verified that there were many differentially enriched immune cell-related pathways between the high-risk group and the low-risk group ( ). Then, 314 DEGs were identified via DESeq2 and edgeR package in the R programme ( ). Further KEGG pathway enrichment analysis indicated that these DEGs were mainly enriched in immune-related pathways ( , Ras signaling pathway, PPAR signaling pathway, and IL-17 signaling pathway). To further study the correlation between FeSig and immune status, we quantified the enrichment scores of diverse immune cell subpopulations, related functions, or pathways with ssGSEA. As shown in , nine immune cell subpopulations (CD8+ T cells, macrophages, mast cells, neutrophils, T helper cells, Tfhs, Th1 cells, Th2 cells, and TILs) were clearly different between the high-risk group and the low-risk group. For immune-related functions, seven pathways (CCR, checkpoint, cytolytic activity, inflammation promotion, parainflammation, T cell coinhibition, and T cell costimulation) had higher scores in the high-risk group, and only the type II IFN response had higher scores in the low-risk group ( ). These findings revealed that FeSig is significantly associated with the regulation of tumor immunity.

Figure 4

The FeSig was obviously closely associated with immune-related pathways. (A–C) GSEA and GSVA showed that FeSig was correlated to immune cell related pathways. (D, E) A total of 314 DEGs were identified via DESeq2 and edgR package. (F) The DEGs were mainly enriched in immune-related signaling pathways. FeSig, ferroptosis-related gene signature; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; DEGs, differentially expressed genes.

Figure 5

ssGSEA analysis of FeSig. (A, B) The ssGSEA score of 16 immune cells and 13 immune-related functions between high-risk group and low-risk group. ssGSEA, single sample gene set enrichment analysis. ***P < 0.001; **P < 0.01; *P < 0.05; ns, no significance.

Potential of the FeSig as an Indicator of Response to Anti-PD-1 Therapy

Due to the wide popularization of immune checkpoint inhibitors in ccRCC, we further focused on the correlation between FeSig and immune checkpoint pathways. As shown in , we found that PD-1, CTLA4, LAG3, and TIGIT were upregulated in the high-risk group and TIM-3 was downregulated in the high-risk group. Moreover, correlation analysis indicated that the riskScore was significantly positively correlated with PD-1, CTLA4, LAG3, and TIGIT ( ). These findings suggested that FeSig may be a biomarker of the response to immunotherapy. Therefore, the predictive value of FeSig was tested in an immunotherapy dataset of ccRCC. Survival analysis verified that there was a significant difference (log-rank p = 0.038) between the high- and low-risk groups of patients with anti-PD-1 monotherapy. Patients in the high-risk group had longer OS than those in the low-risk group ( ). The ROC curve of ICI response ( ) revealed that the diagnostic value of FeSig (AUC = 0.603) was better than that of PD-1 (AUC = 0.544).

Figure 6

FeSig was a potential indicator of response to anti-PD-1 therapy. (A) The expression level of immune checkpoints between high-risk group and low-risk group. (B–H) The correlation between immune checkpoints and risk score. (I) High risk group of patients with anti-PD-1 monotherapy had better OS. (J) ROC curve of ICI response showed that FeSig had good diagnostic value. OS, overall survival; ROC, receiver operating characteristic; ICI, immune checkpoint inhibitor. *P < 0.05; ns, no significance.

Identification of the TAZ/WNT10B Axis as an Immune Checkpoint Regulatory Pathway in ccRCC

The WGCNA algorithm was performed to find further downstream targets of four-gene signature (FeSig). As shown in , β = 5 was selected as the soft threshold power in this study. The DEGs between high-risk group and low-risk group were divided into four modules according to the expression level of each gene ( ). Combined with the clinical information, we found that the MEturquoise module was significantly positively associated with grade (R = 0.23, p = 2e-07), TNM stage (R = 0.29, p = 3e-11), risk (R = 0.66, p = 2e-66), and ImmuneScore (R = 0.2, p = 6e-06) ( ). According to cor.geneTraitSignificance >0.2 and cor.geneModuleMembership >0.6, CPNE7, WNT10B, ADAMTS14, and RUFY4 were regarded as downstream hub genes of FesSig ( , ). Moreover, we analyzed the correlation among four-gene signature (FeSig) (BID, TAZ, MT1G, and SLC7A11), downstream hub genes (CPNE7, WNT10B, ADAMTS14, and RUFY4), and immune checkpoints (PD-1, CTLA4, LAG3, and TIGIT) to further discover potential molecular mechanisms of tumor immunity. According to the correlation value, BID, TAZ, CPNE7, WNT10B, and RUFY4 were selected for subsequent analysis ( ). Survival analyses of OS and DFS proved that WNT10B had better prognostic prediction performance than CPNE7 and RUFY4 ( ). Similarly, TAZ had better prognostic value than BID ( ). Therefore, TAZ and WNT10B were selected for further study. In addition, the expression of TAZ was clearly positively correlated with WNT10B ( , R = 0.66, p = 2e-67). Interestingly, GSEA verified that the low TAZ group was enriched in WNT signaling pathway ( ). These findings indicated that TAZ might regulate tumor immunity through mediating WNT10B (WNT signaling pathway) in ccRCC.

Table 1

Hub genes were identified via WGCNA.

Gene	cor.geneModuleMembership	cor.geneTraitSignificance
	Turquoise	Grade	TNM	Risk	ImmuneScore
CPNE7	0.63	0.29	0.27	0.49	0.41
WNT10B	0.81	0.21	0.21	0.55	0.31
ADAMTS14	0.66	0.29	0.26	0.46	0.22
RUFY4	0.75	0.20	0.24	0.49	0.40

WGCNA, weighted gene co-expression network analysis.

Figure 8

TAZ/WNT10B axis was identified as a potential immune checkpoint regulatory pathway. (A–F) WNT10B had better prognostic value than CPNE7 and RUFY4 according to survival analysis of OS and DFS. (G, H) Survival analysis revealed that TAZ had better prognostic value than BID. (I) TAZ was closely positively correlated to WNT10B (R = 0.66). (J, K) Low TAZ group was enriched in WNT signaling pathways. OS, overall survival; DFS, disease-free survival.

TAZ/WNT10B Was Closely Correlated With TNM/Grade Stage and UpRegulated in ccRCC Tissue

We further analyzed the correlation between the expression level of TAZ/WNT10B and TNM/Grade stage and verified their mRNA expression level in ccRCC tissues. As shown in , TAZ and WNT10B were both elevated and positively correlated with TNM/Grade stage in ccRCC. Moreover, qRT-PCR assay also indicated that TAZ and WNT10B were upregulated in ccRCC tissues ( ).

Figure 9

TAZ/WNT10B was closely correlated with TNM/Grade stage and up-regulated in ccRCC tissue. (A) In TCGA-KIRC dataset, TAZ was elevated in ccRCC and positively correlated with TNM/Grade stage. (B) In TCGA-KIRC dataset, WNT10B was elevated in ccRCC and positively correlated with TNM/Grade stage. (C) The mRNA expression level of TAZ/WNT10B was upregulated in ccRCC tissues. ccRCC, clear cell renal cell carcinoma. ***P < 0.001; **P < 0.01; *P < 0.05.

Discussion

Recently, immunotherapy has emerged as a promising approach for cancer treatment, which mainly includes non-specific immunostimulation, immune checkpoint inhibitors (ICIs), tumor vaccines, and adoptive cellular immunotherapy (28). In metastatic RCC (mRCC), immune checkpoint inhibitors have changed the treatment paradigm because most patients with newly diagnosed mRCC are now treated with these medicines. It has been reported that immune checkpoint inhibitors have provided significant clinical benefit for mRCC patients in multiple clinical trials (29). However, there are still some patients without good effects due to adverse reactions and drug resistance to ICIs (30). Therefore, it is necessary to study the molecular mechanisms of drug resistance and reduce the adverse reactions to ICIs. At present, some studies have reported that the induction of ferroptosis enhances the effect of ICIs (16). Interestingly, many ferroptosis-related genes are abnormally expressed in ccRCC and might be potential therapeutic targets. However, the correlation between ferroptosis and immune checkpoints in ccRCC is still unclear.

In this study, a ferroptosis-related risk model (FeSig) was constructed through Cox regression analysis. We found that FeSig had good diagnostic/prognostic value and was closely associated with immune checkpoints. Moreover, patients with high-risk scores had better OS than those with low-risk scores after PD-1 inhibitor treatment, which indicated that FeSig was a potential prognostic biomarker for immunotherapy. Then, the TAZ/WNT10B axis was identified as a potential regulatory pathway of immune checkpoints via WGCNA and correlation analysis.

TAZ (Tafazzin) is a transcriptional regulator and plays a vital role in tumorigenesis and tumor progression of most solid tumors (31). TAZ stimulates cell proliferation by regulating DNA duplication and mitosis (32). It has been reported that TAZ maintains plasticity in cell–ECM adhesion and favors cytoskeletal remodeling to promote tumor metastasis (33). TAZ is also regarded as an important regulator of ferroptosis. Two research groups recently reported that TAZ promoted tumor cell ferroptosis by regulating members of the NOX family (34, 35). Furthermore, many studies have verified that TAZ is closely involved in tumor immunity and the microenvironment. TAZ mediates the expression of tumor-secreted factors to drive the differentiation and recruitment of immune suppressive cells, such as tumor-associated macrophages (TAMs) (36) and regulatory T cells (Tregs) (37). In addition, TAZ promotes immune evasion of tumor cells by regulating the expression of immune checkpoints. Feng et al. indicated that tumor cell-derived lactate activated the TAZ/PD-L1 axis to enhance tumor evasion from the immune response (38). Similarly, a recent study also reported that the TAZ/YAP/TEAD signaling pathway increased PD-L1 promoter activity and induced immune evasion of tumors (39). In our study, we found that the expression of TAZ was highly positively associated with PD-1, which suggested that TAZ was a potential therapeutic target of immunotherapy in ccRCC.

To study the specific molecular mechanisms of TAZ, WGCNA and GSEA were performed and indicated that WNT10B is a potential downstream target of TAZ. WNT10B (Wnt Family Member 10B) is a regulator encoding secreted proteins and activating the Wnt signaling cascade (40). It has been reported that the Wnt signaling pathway is closely involved in regulating immune checkpoints. Notably, the expression of PD-L1 has been proven to be regulated by MYC, which is a well-documented target of the Wnt signaling pathway (41). Moreover, inhibiting the Wnt/β-catenin axis can promote antitumor immunity by suppressing PD-L1 expression (42). In this study, we also found that TAZ was significantly positively correlated with WNT10B and that high WNT10B predicted poor OS/DFS in ccRCC. These findings suggested that the TAZ/WNT10B axis might regulate tumor immunity by activating the WNT signaling pathway.

In conclusion, the TAZ/WNT10B axis (ferroptosis-related pathway) is regarded as a tumor immune-related regulatory pathway via integrated bioinformatics analysis. Further analysis revealed that immune checkpoints are potential targets of the TAZ/WNT10B pathway. Therefore, the TAZ/WNT10B axis is expected to be a novel therapeutic target of immunotherapy in ccRCC. However, the specific mechanisms still require further research.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: TCGA-KIRC, https://www.cancer.gov/tcga and ICGC-RECA-EU, https://dcc.icgc.org https://science.sciencemag.org/content/suppl/2018/01/03/science.aan5951.DC1.

Ethics Statement

The studies involving human participants were reviewed and approved by Human Research Ethics Committee of Huazhong University of Science and Technology (HUST). The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

XZ designed this study. TX and SG performed data collection and analysis. TX, SG, HR, and JL performed the majority of the experiments. TX and SG wrote the manuscript and contributed to preparing and making figures and tables. YL, DL, and JT collected the clinical samples and managed the clinical data. JT and JS reviewed the relevant literature. HY and XZ provided conceptual advice and critically reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Key Research and Development Plan in China (grant no. 2017YFB1303100), the National Natural Science Foundation of China (grant nos. 81672524, 81672528, and 81874090), the Hubei Provincial Natural Science Foundation of China (grant no. 2018CFA038), the Independent Innovation Foundation of Huazhong University of Science and Technology (grant no. 118530309), and the Clinical Research Physician Program of Tongji Medical College, Huazhong University of Science and Technology (grant no. 5001530015).

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.