Inflammation-Related Gene Signature: An Individualized Risk Prediction Model for Kidney Renal Clear Cell Carcinoma

Background There is much evidence that confirms the inextricable link between inflammation and malignancy. Inflammation-related regulators were involved in the progression of kidney renal clear cell carcinoma (KIRC). However, the predictive role of single gene biomarkers is inadequate, and more accurate prognostic models are necessary. We undertook the current research to construct a robust inflammation-related gene signature that could stratify patients with KIRC. Methods The transcriptome sequencing data along with clinicopathologic information of KIRC were obtained from TCGA. A list of inflammation-related genes was acquired from the Molecular Signatures Database. Using the RNA-seq and survival time data from the TCGA training cohort, an inflammation-related gene signature was built using bioinformatic methods, and its performance in predicting patient prognosis was assessed by Kaplan–Meier and ROC curve analyses. Furthermore, we explored the association of risk score with immune score, stromal score, tumor immune-infiltrating cells (TIICs), immunosuppressive molecules, m6A regulators, and autophagy-related biomarkers. Results Herein, nine inflammation-related hub genes (ROS1, PLAUR, ACVR2A, KLF6, GABBR1, APLNR, SPHK1, PDPN, and ADORA2B) were determined and used to build a predictive model. All sets, including training set, four testing sets, and the entire TCGA group, were divided into two groups (low and high risk), and Kaplan–Meier curves all showed an adverse prognosis for patients in the high-risk group. ESTIMATE algorithm revealed a higher immune score in the high-risk subgroup. CIBERSORT algorithm illustrated that the high-risk group showed higher-level immune infiltrates. Furthermore, LAG3, TIGIT, and CTLA4 were overexpressed in the high-risk subgroup and positively associated with risk scores. Moreover, except for METTL3 and ALKBH5, the other m6A regulators decreased in the high-risk subgroup. Conclusions In conclusion, a novel inflammation-related gene signature comprehensively constructed in the current study may help stratify patients with KIRC.


Introduction
Kidney renal clear cell carcinoma (KIRC) is the most lethal urological tumor and its incidence and mortality are increasing yearly [1]. Radical surgery is the preferred treatment of limited renal clear cell carcinoma. en, 20-40% of patients in the early stages eventually develop metastatic KIRC. Moreover, approximately 30% of patients with renal clear cell carcinoma have a metastasis initial diagnosis due to insidious onset [2]. Unlike other advanced malignancies, advanced renal clear cell carcinoma is resistant to conventional radiotherapy, and although the advent of targeted drugs such as tyrosine kinase and mTOR pathway inhibitors has enhanced the long-term survivals for several patients, the clinical outcome for most patients remains poor due to the presence of toxic side effects and the emergence of drug resistance [3,4]. e link between cancer and inflammation has been explored extensively since it was discovered in the 19th century. Several lines of evidence suggest that tumors usually occur in the site of chronic inflammation and inflammatory cells exist in the biopsy of tumor [5]. Researchers found that inflammation mediators and cellular effects are essential components of the local tumor environment [6]. In several types of cancer, inflammation exists prior to the development of malignant changes. In contrast, carcinogenic changes in other types of cancer can induce an inflammatory microenvironment and promote tumor progress [7]. Whatever its origin, the inflammation in the tumor microenvironment has many tumorigenesis effects. It not only accelerates tumor progression by promoting the proliferation, angiogenesis, and metastasis, but also disrupts adaptive immune responses and makes tumor cells tolerant to hormones and chemotherapy drugs.
is cancer-related inflammatory molecular pathway is now being uncovered [8]. Balkwill et al. [9] have revealed that the invasion ability of neoplastic cells is increased in the presence of inflammatory cytokines. Tan et al. [10] have shown that inflammationrelated genes might serve as important prognostic biomarkers for assessing recurrence risk (GADD45G) and death (CARD9, CIITA, and NCF2) in patients with KIRC. At present, some therapeutic drugs for inflammatory cytokines are being developed and tested in clinical practice [11], suggesting that targeting inflammation-related genes is a promising cancer therapy.
As mentioned above, targeting inflammation-related biomarkers may be a promising novel choice for tumor treatment. A large number of inflammation-related regulators are associated with the KIRC progression; however, cancer is a disease caused by the combined involvement of multiple genes and pathways. Given the limitations of a single biomarker, we screened multiple inflammation-related genes for prognostic relevance and constructed a gene signature for risk stratification and prognostic assessment of patients. Herein, we aim to develop an inflammation-related lncRNA model to predict the survival outcomes of patient with KIRC. We used the TCGA database to develop and validated the individualized prognostic signature for KIRC based on inflammation-related genes. Combined with the inflammation-related genes with clinical variables, we construct a comprehensive gene model that could assess the prognosis of patients with KIRC.

Data Collection.
RNA-Seq gene expression data for KIRC was downloaded from the TCGA database (https:// portal.gdc.cancer.gov/), called TCGA-KIRC. e reads per map per million base pairs (FPKM) counts and fragment counts per thousand transcripts were downloaded for further analysis. We finally obtained RNA sequencing data from 530 patients with complete clinical information and their clinicopathological data.

Identification of Differentially Expressed Inflammation-Related Genes (DE-IFRGs).
A comprehensive list of inflammation-related genes (IFRGs) was retrieved from the hallmark gene sets from the Molecular Signatures Database v7.4 (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp), which consists of 200 IFRGs. e "limma" R package and the Wilcoxon test method were used to identify the DE-IFRGs with an adjusted P < 0.05 between KIRC and adjacent normal renal tissues. e "pheatmap" R package was employed to visualize the degree range of differences in the TCGA-KIRC datasets.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and
Genomes (KEGG). To reveal the potential biological functions and underlying action mechanisms of DE-IFRGs, we conducted the GO and KEGG analyses applying the "clusterProfiler" R package [12]. Functional enrichment items were considered as "functional" when the false discovery rate (FDR) <0.05.
en, we further reduce the amount of genes using the LASSO regression analysis to prevent overfitting. Finally, multivariate assays were conducted to identify the hub IFRGs and build a prognostic signature. We then calculate the risk score for each KIRC patient using the following formula: exp gene 1 * β gene 1 + exp gene 2 * β gene 2 + exp gene 3 * β gene 3 + . . . exp gene n * β gene n. Furthermore, patients in all sets as well as the entire TCGA set were classified into low-and high-risk subgroups according to the median risk score of the training set. en, survival assays were conducted. ROC assays were utilized to measure the predictive capability of the prognostic model.

Evaluation of the Risk Signature. Uni-and multivariate
Cox regression analyses were conducted to select the independent prognostic factors. Besides, the associations between risk scores and clinical features of patients were studied. en, we construct a nomogram consisting of independent prognostic factors to predict the OS of KIRC patients. Calibration curve was employed to compare the differences between predicted OS and actual OS. In addition, we compared the differences in the ability of risk model as well as clinicopathological variables to assess patient prognosis.
2.6. Functional Enrichment Analysis. Differentially expressed genes (DEGs) between the high-and low-risk subgroups were identified using the "limma" R package. Genes with |log2FC| ≥ 1, FDR <0.05 were considered differentially expressed. en, GO and KEGG assays based on these DEGs were carried out applying the "clusterProfiler" R package [12].

Evaluation of the Tumor Microenvironment (TME) and Tumor Infiltrated Immune Cells (TIICs).
e ESTIMATE algorithm was used to evaluate scores representative of the relative proportion of immune and stromal cells. Furthermore, we further compared the difference of immune and stromal scores between high-and low-risk subgroups by the Wilcoxon test. Additionally, to analyze the relationships between risk score and TIICs, the content of TIICs was calculated using the CIBERSORT algorithm (http:// cibersort.stanford.edu/).

Association of Risk Score with Immunosuppressive Molecule, m6A Regulators, and Autophagy-Related Biomarkers.
Considering immune checkpoint inhibitors (ICIs) were clinically employed to treat KIRC, we evaluated the association between risk score with ICI-related regulators. m6A regulators and autophagy-related biomarkers were closely related to cancer progression; we thus evaluated the correlation between risk score and m6A regulators as well as autophagy-related biomarkers.
2.9. Statistical Analysis. All statistical analyses were carried out using R (version 3.6.1). Univariate, LASSO, and multivariate assays were used to select the prognostic genes and develop a gene signature. e Kaplan-Meier analysis was applied to show the survival difference. ROC assays were applied to estimate the predictive performance of the risk model. e independent prognostic factors were determined applying multivariate assays. Wilcoxon's test and Pearson's correlation methods were utilized to evaluate the association of risk score with TME, TIICs, ICI-related regulators, m6A regulators, and autophagy-related biomarkers. P-value <0.05 was considered statistically significant.

Data Preparation.
e detailed workflow flowchart of this study is listed in Figure 1. e transcriptome profiles and clinical information of 530 patients with KIRC were publicly downloaded from the TCGA database. We then randomly divided all patients into the training set (n � 320), testing-1 set (n � 53), testing-2 set (n � 52), testing-3 (n � 52), and testing-4 set (n � 53). Data from the training set was used to choose prognosis-related hub IFRGs and construct a risk signature. Simultaneously, data from testing-1, testing-2, testing-3, and testing-4 sets as well as the entire group was utilized to demonstrate the capability of the risk score.

Identification of DE-IFRGs.
e "limma" R package was employed to screen the differentially expressed DE-IFRGs between KIRC samples and normal renal specimens. Herein, 177 dysregulated genes were identified, of which 46 were downregulated, and 131 were upregulated (Figure 2(a)). Figure 2(b) shows the top ten up-and downregulated IFRGs in KIRC. Additionally, we calculate the Pearson coefficients DE-IFRGs, and Figure 2(c) showed a strongly correlated DE-IFRGs association map (cor > 0.8 and P < 0.05), of which the strongest correlations were found between CXCL11 and CXCL10, LTA and LCK, and MSR1 and C3AR1 (Figures 2(d)-2(f )).

Functional Enrichment Analysis of DE-IFRGs.
Functional enrichment analysis of these DE-IFRGs was conducted using the "clusterProfiler" R package. As revealed in Figure 3(a), the significantly enriched BP terms were response to molecule of bacterial origin, response to lipopolysaccharide, and positive regulation of cytokine production; in terms of CC, DE-IFRGs were mainly involved in positive regulation of cytokine production, secretory granule membrane, and membrane raft; as for MF, DE-IFRGs were mainly involved in receptor-ligand activity, cytokine receptor binding, and cytokine activity. Figure 3(b) showed the three significantly enriched GO terms and relevant DE-IFRGs involved in them. Additionally, the top 10 KEGG pathways were TNF signaling, lipid and atherosclerosis, JAK-STAT signaling, chemokine signaling pathway, Influenza A, Toll-like receptor signaling pathway, and inflammatory bowel disease (Figure 3(c)). Figure 3(d) displays the three significantly enriched signaling pathways and related DE-IFRGs involved in these pathways.

Construction and Validation of a Risk Signature Based on
Prognosis-Related IFRGs. Using univariate Cox regression analysis, 20 prognosis-related IFRGs were identified (P < 0.001) ( Table 1). Subsequently, the least absolute Lasso regression analysis was employed to prevent the overfitting and determine the most important prognosis-related IFRGs in KIRC (Figures 4(a) and 4(b)). en, stepwise multivariate assays were applied to build a gene signature. Eventually, nine hub IFRGs (ROS1, PLAUR, ACVR2A, KLF6, GABBR1, APLNR, SPHK1, PDPN, and ADORA2B) were used to construct the gene signature ( Figure 4(c)). Based on regression coefficients (Table 2), we calculated the risk score for each patient using the following formula: risk score � (1.069 * ROS1) + (0.339 * PLAUR) + (−0.720 * ACVR2A) + (−0.198 * KLF6) + (0.600 * GABBR1) + (−0.164 * APLNR) + (−0.386 * SPHK1) + (0.183 * PDPN) + (0.472 * ADORA2B). As exhibited in Figures 5(a) and 5(b), ROS1, PLAUR, GABBR1, SPHK1, and PDPN were overexpressed in the high-risk subgroup, whereas ACVR2A, KLF6, and APLNR were distinctly decreased in the high-risk subgroup. However, no difference was found in ADORA2B. Moreover, our group observed that overexpression of ROS1 and PLAUR indicated worse overall survival (Figures 5(c) and 5(d)). e downregulation of ACVR2A and KLF6 predicted a poor prognosis of patients (Figures 5(e) and 5(f )). Increased expression of GABBR1 was associated with a shorter OS ( Figure 5(g)). Low APLNR expression predicted a shorter OS ( Figure 5(h)). Increased expression of SPHK1 and PDPN suggested worse prognosis (Figures 5(i) and 5(j)). Also, no difference was found in ADORA2B ( Figure 5(k)). Furthermore, using the cBio-Portal database, we explored the genetic mutations of 9 hub IFRGs, and results were shown in Figure 5(l). Subsequently, 320 patients in the training set were stratified into the lowand high-risk subgroups based on the median risk score value. Kaplan-Meier curves showed that high-risk patients showed poorer OS by comparison with low-risk patients (P < 0.001) (Figure 6(a)). ROC assays were utilized to evaluate the prognostic performance of the gene signature, and results showed that the area under the ROC curve for 1-year, 3-year, and 5-year OS was 0.766, 0.721, and 0.751 ( Figure 6(b)). e survival status and the expressions of 9-IFRGs in the training cohort were presented in −. To verify the predictive performance of the gene model, patients in the testing-1 cohort, testing-2 cohort, testing-3 cohort, testing-4 cohort, and the entire group were classified as high-and low-risk subgroups. e Kaplan-Meier survival curve showed a significantly good OS in the lowrisk group (Figures 7(a)-7(e)). e AUC of the gene signature in the testing-1 cohort for 1-year, 3-year, and 5-year OS is also shown in Figures 7(f )-7(j).

Independent Prognostic Analysis, Correlation of Risk Score with Clinical Features, and Construction of a Nomogram.
By coupling with the risk model and clinicopathological features, we identified the risk score (HR � 1.023, P < 0.001) as a factor of overall survival for KIRC using uni-and multivariate Cox regression analyses (Figures 8(a) and 8(b)). Besides, we showed that elevated risk score was notably correlated with higher histological grade (P < 0.05, Figure 8(c)), advanced clinical stage (P < 0.05, Figure 8(d)), and T stage (P < 0.05, Figure 8(e)), suggesting that risk score was positively correlated with tumor progression. Moreover, we used the independent prognostic factors to establish a prognostic nomogram (Figure 8         of DEGs, respectively. As illustrated in Figure 9(c), concerning biological processes, DEGs were significantly enriched in the modulation of negative regulation of hydrolase activity; with regard to cellular components, DEGs were significantly involved in the collagen-containing extracellular matrix, presynapse, and synaptic membrane; in point of molecular functions, DEGs were noticeably involved in receptor-ligand activity, passive transmembrane     transporter activity, and channel activity. DEGs were mainly enriched in phototransduction, linoleic acid metabolism, cholesterol metabolism, arachidonic acid metabolism, IL-17 signaling pathway, and protein digestion and absorption (Figure 9(d)).

3.7.
Association of Risk Score with TME. Immune and stromal cells are crucial constituents of the immune microenvironment. In this current study, the contributions of stromal and immune cells to KIRC were estimated by the ESTIMATE algorithm. e results signified that immune score was crucially higher in the high-risk group (Figure 10(a)); however, no difference was found for the stromal score (Figure 10(b)). Additionally, we applied the CIBERSORT algorithm to compare the differences in each type of immune infiltrating cells. Figure 10(c) showed the proportion of 21 immune cells in each sample. Figure 10(d) illustrates the correlations between infiltrated immune cells in the tumor. Figure 10(e) shows the heat map of the 21 immune cell proportions. Moreover, Figure 10(f ) shows the relationship between risk score with different immune cells, and we found that the high-risk group showed higher-level immune infiltrates of M0 macrophages, regulatory T cells (Tregs), follicular helper T cells, plasma cells, and memory B cells.

Association of Risk Score with Immunosuppressive Molecules, m6A Regulators, and Autophagy-Related Biomarkers.
en, we estimated the association between immunosuppressive molecules and risk score. Figure 11(a) shows the heat map of common immunosuppressive molecules in high-and low-risk subgroups. Furthermore, as illustrated in Figure 11(b), patients with high-risk score expressed higher levels of LAG-3, ICOS, CTLA4, PDCD1, CD27, and TIGIT, whereas HAVCR2 was overexpressed in patients with the low-risk score. Correlation analysis confirmed that LAG-3 (cor � 0.15, Figure 11(c)), TIGIT (cor � 0.11, Figure 11(d)), and CTLA4 (cor � 0.21, Figure 11(e)) were positively associated with the risk score, whereas no difference was found for ICOS, PDCD1, CD27, and HAVCR2 (Figures 11(f )-11(i)). Together, these results indicate that LAG-3, TIGIT, and CTLA4 were positively associated with the risk score. Recent evidence indicated the vital role of m6A mRNA methylation in reducing the antitumor response of CD8 + T cells and promoting anti-PD-1 drug resistance. We thus assess the relationship between risk score and m6A regulators. Figure 12(a) shows the heat map of common m6A regulators in high-and low-risk subgroups. Additionally, we discovered that most of the m6A regulators were significantly decreased in the high-risk subgroup except for METTL3 (Figure 12(b)).
e results indicate that high-risk subgroup patients may be more suitable for immunotherapy with emerging checkpoint inhibitors. Growing researches have revealed a key role for autophagic pathways and proteins in immunity and inflammation. We thus explore the association of autophagy-related genes with risk score, and we found that several autophagy-related genes have a significant link with risk score (Figure 13(a)), and the top three relevant autophagy-related genes are DKK1 (Figure 13(b)), SNAI2 (Figure 13(c)), and AREG (Figure 13(d)).

Discussion
In this work, we constructed an inflammation-related gene feature and evaluated its predictive capability in predicting OS of KIRC patients. en, we studied the potential functions and signaling pathways closely related to risk score and further explored the association between risk score with immune microenvironment, immunosuppressive molecules, m6A regulators, and autophagy-related biomarkers.
Here, nine hub IFRGs (ROS1, PLAUR, ACVR2A, KLF6, GABBR1, APLNR, SPHK1, PDPN, and ADORA2B) were selected by bioinformatics and used to construct 9-IFRG risk signature successfully. Afterwards, we found that gene signature performed well in the training set, testing-1 set, testing-2 set, testing-3 set, testing-4, and the entire TCGA group. Specifically, the higher the risk score of patients is, the worse the overall survival rate is. ROC curve also confirms the robust predictive performance of the risk model. Additionally, by combining the risk model with the clinicopathological features of patients, we found that the 9-IFRG gene model can independently predict the OS of patients with KIRC. Further investigation indicated that the nomogram performed well at predicting 1-, 3-, and 5-year OS in KIRC patients. Furthermore, we found that the risk score was significantly associated with cancer progression in KIRC patients. Moreover, compared to other clinical variables, the risk score had the highest predictive performance of prognosis. To sum up, we constructed a powerful 9-IFRG risk signature and an effective nomogram for KIRC risk stratification and overall survival prediction.
Of the nine hub IFRGs (ROS1, PLAUR, ACVR2A, KLF6, GABBR1, APLNR, SPHK1, PDPN, and ADORA2B) we identified, some are associated with cancer progression. e protooncogene ROS1 encodes a tyrosine kinase receptor that has an essential physiological role in humans. Studies have shown that somatic chromosomal fusions involving ROS1 generate chimerical tumor proteins that can cause various     cancers [13]. In inflammatory myofibroblastic tumors, ROS1 expression predicts ROS1 gene rearrangement [14]. PLAUR, also known as u-PAR, is an essential molecule in modulating cell surface fibrinogen activation and plays a vital role in many healthy and pathological processes [15]. Abnormal PLAUR disorders played a key role in the progression and metastasis of human colon cancer [16]. Moreover, PLAUR impacted colorectal liver metastases by influencing the protein hydrolytic activity and inflammation of the tumor microenvironment in colorectal cancer. Consequently, the colorectal liver metastases [17] ACVR2A is a ligand for activin A protein and is closely associated with polyarthrosis syndrome, protointestinal embryogenesis, and spermatogenesis [18]. Emerging evidence indicated that ACVR2A is involved in many cancer-related signaling pathways, such as the PEDF-induced signaling, the TFG-β signaling, or signaling pathways regulating stem cell pluripotency [19]. KLF6 is a transcription factor of the zinc finger family and modulates lipid homeostasis in KIRC [20]. Additionally, KLF6 had been found to promote the expression and function of proinflammatory genes by inhibiting miR-223 expression in macrophages [21]. GABBR1, also known as GABABR1, is a 7-transmembrane receptor. In colorectal cancer, decreased GABBR1 fosters the proliferation and invasion; overexpression of GABBR1 has the opposite [22]. APLNR is also a seven-transmembrane G protein-coupled receptor that is universally present in diverse tissues. In osteosarcoma, elevated APLNR expression promotes proliferation and invasion [23]. SPHK1 is a biologically active metabolite of sphingosine that is involved in various tumor progression by enhancing cell proliferation and motility. Currently, drugs targeting SPHK1 are now being progressively validated in clinical trials [24]. Type I integral membrane glycoprotein encoded by PDPN is widely distributed in human tissues. In breast tumor-infiltrating immune cells, PDPN was found highly expressed in tumorassociated macrophages (TAMs), and the latter spurs local stromal remodeling and promotes vascular growth and lymphatic infiltration [25]. ADORA2B is a member of the G protein-coupled receptor superfamily and encodes an adenosine receptor. A recent report indicates that hypoxiainducible factor 1-dependent expression of ADORA2B facilitates breast cancer stem cell enrichment [26]. e above reports confirmed the role of 9 hub IFRGs in carcinogenesis. However, whether ROS1, PLAUR, ACVR2A, KLF6, GABBR1, APLNR, SPHK1, PDPN, and ADORA2B affect the clinical outcome of KIRC patients via modulating the process of inflammation requires to be further elaborated, and there are few relevant studies.
To elucidate the functional roles associated with the risk score, the DEGs between the high-risk and low-risk subgroups were identified and used to perform functional enrichment analysis. Intriguingly, we noticed that DEGs are involved in several tumor-related signaling pathways. ese signaling pathways are all in connection with the regulation of tumor immunity.
rough the interaction between chemokines or cytokines and their receptors, different subsets of immune cells are recruited into the tumor microenvironment, causing these populations having a differential impact on tumor progression and treatment outcome [27]. In gastric cancer, elevated intratumoral mast cells resulted in immune suppression via modulating TNFα-PD-L1 pathway [28]. e JAK-STAT signaling pathway is involved in tumor cell recognition and tumor-driven immune escape and plays a role in almost all immune regulatory processes [29]. Toll-like receptor signaling pathway is a classical immune signaling pathway that plays an irreplaceable role in modulating tumor immunity and cancer progression [30]. In addition, we found that the high-risk group had a higher immune score. With regard to immune infiltrating cells, we found that high-risk group showed higher level immune infiltrates. Among them, regulatory T cells (Tregs) play crucial roles in keeping self-tolerance and immune homeostasis. However, in some cases, they promote tumor progression by inhibiting the effective antitumor response [31]. e low-risk group showed higher level immune infiltrates. Among them, M1 macrophage types are thought to be key factors in antitumorigenesis, production of proinflammatory cytokines, and promotion of T-cell immunity. Furthermore, the study suggested that LAG-3, CTLA4, and TIGIT were highly expressed in the high-risk subgroup and also positively associated with risk score, indicating that the high-risk group is in a more immunosuppressed state by comparison with the low-risk group, but also means that patients of high-risk subgroup may benefit more from immune checkpoint inhibitors. N6methyladenosine (m6A) RNA methylation plays a crucial role in the tumor immune microenvironment cancer development. A recent study indicates that downregulated m6A-related genes predict unfavorable outcomes in gastric cancer [32]. We assess the association of m6A regulators with risk score, and we found that most of the m6A regulators were significantly decreased in the high-risk subgroup. Autophagy is an essential homeostatic process by which cells decompose their components. Recent studies have uncovered a key role for autophagic pathways and proteins in immunity and inflammation. We thus evaluate the association of autophagy-related genes and the risk score, and results indicate that many autophagy-related genes were significantly correlated with risk scores, particularly the DKK1, SNAI2, and AREG.

Conclusion
Collectively, our study constructs and validates a robust 9-IFRG risk signature, which may be to the advantage of risk classification and prognosis prediction in KIRC patients. However, there are still some restrictions that should not be overlooked. Our results are mainly derived from bioinformatic analysis; clinical samples and cellular experiments are required to prove our findings; in addition, our analysis discovered that inflammation-related genes might influence renal clear cell carcinoma progression through several mechanisms; nevertheless, further in vivo and in vitro experiments are needed to explore the exact biological roles.
Data Availability e datasets used and/or analyzed during the present study are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.

Authors' Contributions
Ze Zhang and Yan-Yan Wei contributed equally to this work.