Immune ULBP1 is Elevated in Colon Adenocarcinoma and Predicts Prognosis

Background: Colon adenocarcinoma (COAD) is still the main cause of cancer deaths worldwide. Although immunotherapy has made progress in recent years, there is still a need to improve diagnosis, prognosis, and treatment tools. UL-16 binding protein 1 (ULBP1) is a ligand that activates the receptor natural killer cell group 2 receptor D (NKG2D) and plays an important immunomodulatory role. We aimed to investigate the clinical significance of ULBP1 in COAD. Methods: We obtained the relevant data from The Cancer Genome Atlas (TCGA). A total of 438 patients with COAD were included in this study, with a mean age of 67.1 ± 13.03 years old, of which 234 (53.42%) were male. The diagnostic value of COAD tumor tissues and adjacent tissues was analyzed by ROC curve. Univariate and multivariate survival analysis investigated the prognostic value of ULBP1 gene, and Gene Set Enrichment Analysis (GSEA) curve was performed to analyze the biological process and enriched enrichment pathway of ULBP1 in COAD. Combination survival analysis investigated the combined prognostic effect of prognostic genes. Results: ULBP1 gene had a high diagnostic value in COAD [AUC (TCGA) = 0.959; AUC (Guangxi) = 0.898]. Up-regulated ULBP1 gene of patients with COAD predicted a worse prognosis compared to those patients with down-regulated ULBP1 gene (Adjusted HR = 1.544, 95% CI = 1.020–2.337, p = 0.040). The GSEA showed that ULBP1 was involved in the apoptotic pathway and biological process of T cell mediated cytotoxicity, regulation of natural killer cell activation, and T cell mediated immunity of COAD. The combination survival analysis showed that the combination of high expression of ULBP1, AARS1, and DDIT3 would increase the 2.2-fold death risk of COAD when compared with those of low expression genes. Conclusion: The immune-related ULBP1 gene had diagnostic and prognostic value in COAD. The combination of ULBP1, AARS1, and DDIT3 genes could improve the prognostic prediction performance in COAD.


INTRODUCTION
According to Global Sung et al. (2021), it is estimated that there will be more than 1.9 million new cases of colorectal cancer (CRC) and 935,000 deaths, accounting for about one-tenth of cancer cases and deaths. Overall, the incidence of CRC ranks third, but the mortality rate ranks second (Sung et al., 2021). Colon adenocarcinoma (COAD) is a common malignant tumor of the digestive system, and it is the most frequently diagnosed histological subtype of CRC (Siegel et al., 2020). The patient's clinical symptoms usually manifest as diarrhea, abdominal pain, and bloody stools, which develop in the middle and late stages of the disease. The quality of life of patients is usually very low, and the prognosis of most patients is poor. The occurrence and development of COAD are the results of a variety of mixed factors in vivo and in vitro, which involve a series of molecules and signal pathways . Although substantive diagnosis and treatment strategies such as surgery, neoadjuvant therapy, and targeted therapy are constantly being developed, the recurrence rate of postoperative COAD is still high, and the 5year survival rate of patients with advanced COAD is still very low (The 5-year survival rate after distant metastasis is less than 15%) (Patel et al., 2021). Thus, there is an urgent need to explore new biomarkers and therapeutic targets in clinical practice to improve the survival rate of patients with COAD.
Natural killer cell group 2 receptor D (NKG2D) is an alkalineactivated receptor belonging to the c-type lectin-like family. It is expressed in NK cells, most NKT cells, some γδ T+, and CD8 T+ cells. It is different from other NKG2 receptors, which is not associated with CD94 (Mondelli, 2012) and has nothing to do with CD94 (Mondelli, 2012). The seemingly unchanged activation receptor NKG2D is mixed with a variety of ligands, such as the major histocompatibility complex class I-related chain A and B (MICA and MICB) and the unique long 16 (UL16)binding protein family (ULBPs, ULBP1-6) which are poor (Champsaur and Lanier, 2010). NKG2D ligand expression is usually lacking in healthy tissue but can induce expression under stress, infection, and DNA damage. NKG2D ligand is also widely expressed in a variety of cancer cell lines and primary solid tumors (McGilvray et al., 2009;Champsaur and Lanier, 2010). The upregulation of these ligands may break NK cells from inhibiting the balance of activation (induced selfidentification), with significant biological significance. The interaction of NKG2D is variable between different types of cancers. In the mouse model, the tumor cell line of transfection of RaE1 is rejected by NKG2D-mediated immunization (Diefenbach et al., 2001). The most recent NKG2D knockout mice provide the most convincing evidence for NKG2D to participate in anti-tumor immune responses (Guerra et al., 2008). Many mechanisms have been proposed, cancer can evade NKG2D-mediated immune response. In some systems, the persistent expression of NKG2D ligands can cause NKG2D expression to be lowered (Oppenheim et al., 2005). These results indicate that NKG2D's participation in the anticancer immune response is significantly different between different types of cancer. It is also proposed that tumors may release soluble NKG2D ligands, or secrete immunosuppressive cytokines, such as transforming growth factor-beta to reduce NKG2D expression (Groh et al., 2002;Castriconi et al., 2003;Lee et al., 2004). Notably, NKG2D ligands can be independently expressed in human cell lines and primary tumors, the expression of NKG2D ligands among different tumors in knockout mice is also heterogeneous (McGilvray et al., 2009). The complex interaction between NKG2D and its ligands may involve the natural history and treatment response of various cancers (Mondelli, 2012).
The authors showed that ULBP1, one of the important ligands of NKG2D, is up-regulated in COAD cancer tissues, but is lowexpressed in normal adjacent tissues. Although most previous studies reported that ULBP1 was related to recurrence-free survival, disease-free survival, or overall survival (OS) in different cancers (McGilvray et al., 2009;McGilvray et al., 2010;Mondelli, 2012;Chen et al., 2013;Maccalli et al., 2017), the relationship between ULBP1 and OS in COAD has not been reported yet. Therefore, our study uncovers and investigates the diagnosis, prognosis, and immune mechanism of ULBP1 gene in COAD, which may help make this immune receptor an exceptional candidate for basic and applied cancer research in COAD.

Public Data Collection
We downloaded the COAD-related ULBP1 gene mRNA expression data set and the corresponding patient clinical information parameters from the public cancer database-The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/, obtained on December 10, 2020) (Giordano, 2014;Hutter and Zenklusen, 2018). Based on the TCGA-COAD project data, the differential expression level of ULBP1 in tumor tissues and adjacent normal tissues in pan-cancers was obtained from the TIMER website (https://cistrome.shinyapps.io/timer/, obtained on June 10, 2021) (Liu et al., 2021). The expression level of ULBP1 gene in COAD tumor tissues and adjacent normal tissues was also obtained from the GEPIA website, which integrated the expression levels of normal tissues in the TCGA database and the Genotype-Tissue Expression (GTEx) database. Additionally, the expression level of ULBP1 in different COAD tumor stages was obtained from Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn/index.html; obtained on June 10, 2021) (Tang et al., 2017). We obtained the methylation level and mutation status of ULBP1 gene from UALCAN database (http://ualcan.path.uab.edu/index.html, obtained on June 10, 2021) (Li et al., 2021) and cBio Cancer Genomics Portal (cBioPortal, https://www.cbioportal.org/, obtained on June 12, 2021), respectively (Harbig et al., 2021). The ULBP1-expressed protein in COAD cancer tissues and adjacent normal tissues was obtained from THE HUMAN PROTEIN ATLAS (HPA, https://www.proteinatlas.org/, obtained on June 12, 2021) (Ullah et al., 2021). Finally, we obtained the information on immune infiltration associated with ULBP1 in COAD and the correlation between ULBP1 gene expression level and Genomics of Drug Sensitivity in In the TCGA database, patients were divided into two high-and low-expression groups based on the median cut-off value of ULBP1 gene expression. Univariate and multivariate survival analysis was performed to assess the potential prognostic value of ULBP1 gene expression in patients with COAD. According to the expression level of ULBP1 gene, the COAD expression genome-wide data in the TCGA database was divided into high expression group and low expression group. When the gene expression satisfied |log 2 foldchange|≥1 and p < 0.05, it was considered to be a differential expression gene in this study.
We investigated the co-expression analysis of ULBP1 and COAD-related genes in the TCGA cohort. When the Pearson correlation coefficient ≥0.3 or ≤ −0.3, and the p-value < 0.05, it was considered to be a co-expressed gene with the ULBP1 in COAD. The top 20 co-expressed genes were selected to analyze the prognostic value of genes in COAD. Significant prognostic genes were selected to construct a risk score model based on the prognostic contribution coefficients (β) of different genes. The risk score was generated based on the calculation formula: gene expression of 1*β1+ gene expression of 2*β2+. . .+gene expression of n*βn (Ruan et al., 2020b).

Statistical Analysis
In this study, the unpaired Student's t-test or paired t-test was used to compare the expression levels between two groups. The gene expression level was expressed by using the mean ± standard deviation (SD). If the data did not conform to the normal distribution, the rank sum test was used. Univariate and multivariate cox regression analyses were performed to investigate the prognostic value of genes. The selection of adjustment variables adopted single-factor meaningful clinical parameters, and the TNM stage was used as an adjusted factor for prognostic adjustment to reduce the clinical deviation. All twotailed p < 0.05 were considered statistically significant. The statistical analyses in this study were performed using SPSS 22.0 version and R platform, version 4.0.1.

Baseline Characteristics
In this study, a total of 456 patients with ULBP1 mRNA expression data set in COAD were obtained from the TCGA database, including 480 tumor tissue samples and 41 adjacent normal tissue samples. After removing the duplicate information and the information with a survival time of 0, we obtained a total of 438 tumor sample information and 42 adjacent normal tissue sample information. The mean age of the 438 patients was 67.1 ± 13.03 years old, including 234 males (53.42%) and 204 females (46.58%). Clinical parameter information included age, sex, and TNM stage. The univariate survival analysis of clinical parameters showed that only the TNM stage had a significant prognostic value in COAD (p < 0.001, Table 1).

Investigating the Association Between UL-16 Binding Protein 1 Gene Expression and Immune Infiltration and Drug Sensitivity
Based on the GSCA website, the association between ULBP1 gene expression and immune infiltration suggested that the ULBP1 expression was significantly positively related to the cells of nTreg, iTreg, Neutrophil, Monocyte, Gamma_delta, Exhausted, and CD8_navie. However, the inverse relationship was observed in the cells of NK, NKT, Tfh, Th17, Th2, Tr1, MAIT, Cytotoxic, CD8_T, CD4_T, and B cell. Additionally, a correlational relationship was observed in ULBP1 gene expression and the majority of drug sensitivity ( Figure 1).

Differential Expression Analysis and Diagnostic ROC Curve Analysis
We obtained the ULBP1 gene in COAD tumor tissues and normal tissues from the GEPIA database that matched the information from the GTEx database and found that ULBP1 expression was Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 13 | Article 762514 up-regulated in tumor tissues and also found that the expression level of ULBP1 gene increased with the progression of tumor stages (Figures 2A,B). It was also found that the expression level of ULBP1 in most tumor tissues was higher than that in normal tissues adjacent to cancer ( Figure 3A). There was no significant difference in the methylation level of ULBP1 gene in tumor tissues and adjacent normal tissues in COAD, and the mutation rate of ULBP1 gene in COAD was 0% ( Figures  2C,D). ULBP1 expressed protein was mainly expressed in the cytoplasm ( Figures 2E-G).
Based on the TCGA cohort, we analyzed the expression level and diagnostic value of ULBP1 gene in COAD. The results showed that the expression level of ULBP1 gene in COAD tumor tissues was significantly higher than that in adjacent normal tissues (p < 0.001). Simultaneously, it had a higher diagnostic value in COAD (AUC = 0.898, 95%CI = 0.784-1.000, p < 0.0001) ( Figure 3B).
The validation result based on the Guangxi cohort found that the expression level of ULBP1 gene in COAD tumor tissue was significantly higher than that in adjacent normal tissues (p = 0.0028). The diagnostic ROC curve results showed that ULBP1 has a higher diagnostic value in COAD (AUC = 0.959, 95%CI = 0.942-0.976, p < 0.001) ( Figure 3C).

Survival Analysis of UL-16 Binding Protein 1 Gene in Colon Adenocarcinoma
The results of univariate survival analysis included the TNM stage as an adjusted prognostic factor. After adjustment, the high expression of ULBP1 gene in COAD predicted a worse OS compared to patients with low expression of ULBP1 (Adjusted HR = 1.544, 95%CI = 1.020-2.337, p = 0.040) ( Figure 4A; Table 2).
We divided 438 COAD patients into high-and low-expression groups based on the median cut-off value of ULBP1 expression. At the same time, the COAD genome-wide data was also divided into two groups, and the differential analysis and enrichment pathway analysis of these two groups were carried out. The PHGDH gene in the high expression group was significantly up-regulated, while the down-regulated genes included ITLN1, JCHAIN, DUOXA2, CLCA1, PRAC1, ADH1B, GCG, IGLL5, NXPE4, DUOX2, CHP2, and SI. The enrichment pathway analysis showed that these differential genes might involve in the process of extracellular exosome and immunoglobulin receptor binding ( Figure 5).
Based on the ULBP1 expression levels, we performed the Gene Set Enrichment Analysis (GSEA) to investigate the potential prognosis molecular mechanism of ULBP1 in COAD. The GSEA was performed by the tool of GSEA 4.1.0 version. The internal reference genes of GSEA were obtained from the Molecular Characterization Database (MSIGDB): KEGG pathway: c2.cp.kegg.v7.4.symbols.gmt; GO term: c5.go.v7.4.symbols.gm. In this study, nominal p < 0.05 and false discovery rate (FDR) < 0.25 were considered statistically significant. The results showed that ULBP1 gene might involve in the development of COAD by participating in the apoptosis pathway and the biological process of T cell mediated cytotoxicity, regulation of natural killer cell activation, and T cell mediated immunity. (Figure 6; Supplementary Figure S1).

UL-16 Binding Protein 1 Related Co-expression Analysis and Prognostic Analysis in Colon Adenocarcinoma
Based on all gene expression sequences of the TCGA database, the potential clinical value of ULBP1 gene and ULBP1 related coexpressed genes was investigated. A total of 87 co-expressed genes related to ULBP1 in COAD were mined. Pathway analysis of 87 co-expressed genes showed that co-expressed genes were involved in metabolic pathways of COAD (Figure 7).
The model of risk score we constructed by the formula: Risk score = ULBP1*0.434 + DDIT3*0.442 + AARS *0.459. The higher the gene expression, the higher the risk score, and the higher the patient's risk of death. The time-dependent ROC curve results showed that the 1-, 3-, and 5-year AUCs were 59.2, 56.8, and 57.5, respectively ( Figure 8).

DISCUSSION
It is now widely accepted that tumors develop methods to evade anti-cancer immunity through a process called immunoediting (Dunn et al., 2004). Various evidence from in vivo models indicates that the immune system attacks early-stage tumors. Surviving cancer cells must adapt to avoid the immune system. This process is described as immunoediting, immune sculpting, or cancer immune evasion (Dunn et al., 2004). In recent years, studies on tumor models in vivo strongly indicate that the activated immune receptor NKG2D participates in the anticancer immune response, and it has also attracted much attention as a ligand of the NKG2D receptor (Smyth et al., 2005;Guerra et al., 2008;McGilvray et al., 2009). In humans, primary tumors and tumor cell lines express NKG2D ligands at a high frequency (McGilvray et al., 2009). As an important member of NKG2D ligand, ULBP1 also plays an important role in immune regulation. The expression of ULBP1 was associated with majority of immune cells, including NK cells, most NKT cells, γδ T + and CD8 T+ cells. etc. Interestingly, we can take targeted chemotherapy based on the results of drug susceptibility to the immune-related ULBP1 gene.
After mining through the database, it was found that ULBP1 was highly expressed in the majority of tumors (including COAD), compared with adjacent tumor tissues. However, some tumors expressed the opposite trend, such as lung adenocarcinoma (LAUD) ( Figure 2A). As previously described, it also showed that the expression of ULBP1 in different cancers was different, but the general expression was frequently expressed in cancer tissues. When we analyzed the diagnostic value of ULBP1 gene expression differences in COAD, whether it was the TCGA cohort or the Guangxi validated cohort, we found that ULBP1 had a higher diagnostic value (TCGA cohort:0.959; Guangxi cohort:0.898) in COAD. In other words, we can take advantage of this high expression characteristic in cancer tissues, and the immune-related ULBP1 can better distinguish cancer tissues from normal tissues. Additionally, we also found that with the progress of TMN staging, the expression of ULBP1 showed an upward trend. The expression level of ULBP1 gene was related to the tumor grade and prognosis, and the differential expression level of ULBP1 gene was different in tumor tissues and adjacent normal tissues of patients with different tumor stages. Univariate and multivariate survival analysis results showed that low expression of ULBP1 in patients with COAD had a worse prognosis when compared with those patients with high expression of ULBP1. The differential expression results between the high and low groups of ULBP1 expression indicated that it was related to the binding of immunoglobulins. In addition, the GSEA of ULBP1 gene in COAD suggested that ULBP1 was involved in the occurrence and development of COAD through enrichment of apoptosis pathways, and was related to the immunoregulation of T cells and NK cells. We suspected that upregulated-ULBP1 might participate in the apoptosis process of COAD through its unique immune regulation mechanism. A study by CADOUX et al. also found that ULBP1 was expressed at a higher level in hepatocellular carcinoma (HCC) tumors with lower differentiation and higher grades, but the difference is not significant (Cadoux et al., 2021). Interestingly, a study of 462 primary colorectal tumors by McGilvray et al. investigated the ULBP1-expressed protein in different TNM stages, the result showed the opposite trend was that high expression level of ULBP1 was common in TNM stage I tumors, but gradually decreased in stage II, III, and IV tumors (McGilvray et al., 2009). To understand the difference, we also investigated the expression level of ULBP1 in rectal adenocarcinoma from the data platform (Supplementary Figure S2), the trend was consistent with the description of McGilvray et al. (McGilvray et al., 2009), indicating that the expression of ULBP1 in the colon and rectum was also heterogeneous. Changes in the expression level of ULBP1 are inseparable from tumor differentiation and grade. In other words, ULBP1 is closely related to tumor prognosis. However, previous reports described the potential mechanism and prognosis of ULBP1 expression changed. A study of genome-wide screen to identify novel drivers of ULBP1 expression by Gowen et al. showed that in the multiple stages of ULBP1 biogenesis, independent pathways gradually play a role. The transcription factor ATF4 drives the expression of ULBP1 gene in cancer cells, while the RNA binding protein RBM4 supports the expression of ULBP1 by inhibiting a new alternative splicing subtype of ULBP1 mRNA, and explains its mechanism of activating the body's immune system (Gowen et al., 2015). The study by Chava et al. indicated that DOT1L inhibition could regulate apoptotic and metabolic pathways as well as upregulate the expression of ULBP1 that increased in NK cell-mediated ovarian cancer eradication (Chava et al., 2021). Maccalli et al. showed that patients with melanoma with the negative expression of sULBP-1 were associated with a better prognosis than those patients with positive expression of sULBP-1 (OS: 25.3 months vs. 12.1 months) (Maccalli et al., 2017). On the contrary, a study by CADOUX et al. showed that the high expression of ULBP1 was The interaction between NKG2D and its ligands may play a central role in anti-tumor surveillance. The level of NKG2D ligands may determine the strength of the anti-tumor immune response (Wu et al., 2012). As described above, tumors can lead to tumor re-editing through immune evasion or ligand shedding. Different cancers are heterogeneous, and we should treat different cancers differently in their anticancer immune responses. To directly avoid NKG2D recognition, tumors may secrete TGF-β and/or release soluble NKG2D ligands, thereby down-regulating the expression of NKG2D (Lee et al., 2004). This observation was also observed in NKG2D knockout mice. For example, the incidence of MCA-induced fibrosarcoma was not affected when knocked out, but the incidence of large prostate tumors, when knocked out, was much higher than that of wild-type (Guerra et al., 2008). Butler et al. confirmed that p53 family members play an important role in the upregulation of ULBP1 in head and neck squamous cell carcinoma induced by proteasome inhibitor drugs (Butler et al., 2009). It is well known that the activation of immune response by NKG2D depends on the tissue microenvironment and synergizes/antagonizes the signals induced by other cell receptors and cytokines . A similar description was suggested by CADOUX et al. that the activated NKG2D system led to a strong inflammatory  response, leading to a strong aggressiveness and poor prognosis (Cadoux et al., 2021). These factors vary for different types of cancer. It is also clear that NKG2D ligand can be independently expressed on cancer cells and can be expressed in response to different cancer-related pathways. Such as ULBP1-2, but not ULBP3, is induced by the expression of the BCR/ABL oncogene (McGilvray et al., 2010). The enriched pathways of ULBP1 gene and its co-expressed genes showed that co-expressed genes might participate in the metabolic pathway and Aminoacyl-tRNA biosynthesis of COAD. When we selected the top 20 co-expressed genes and performed prognostic analysis, we found that both ARRS1 and DDIT3 genes have prognostic value in COAD. Notably, the combination of High expression of ULBP1, AARS1, and DDIT3 would increase the 2.2-fold death risk of COAD, when compared with those of low expression genes. AARS1 is a family member of the aminoacyl-tRNA synthetases (AARSs), which is a housekeeping protein widely present in all organisms, it can catalyze the combination of amino acids and tRNA and convert nucleic acid coding information into amino acids, playing an important role in protein synthesis (Zhang et al., 2020). In addition to these translation functions, AARSs are also involved in many other important physiological activities, such as translation and transcription regulation, signal transduction, cell migration, angiogenesis, inflammation, and tumorigenesis (Kim et al., 2011;Datt and Sharma, 2014;Kim et al., 2014). Cancer is a disease of cell disorders, which can be affected by using translation in unexpected ways, using the catalytic function of AARSs in an untranslated environment, or manipulating its regulatory function independent of enzyme activity (Wang and Yang, 2020). If the expression of tRNA exceeds a certain level, it may cause abnormal cell and tissue growth. On the other hand, with the strong demand for protein synthesis by cancer, the classic enzyme action of AARSs is needed to maintain tumor growth (Grewal, 2015). DDIT3 gene, also called CHOP, is an endoplasmic reticulum (ER). This gene encodes a member of CCAAT/enhanced binding protein (C/EBP) family transcription factors (Ron and Habener, 1992). DDIT3, activated by p38 mitogen-related protein kinase, is a major pro-apoptotic transcription factor induced by ER stress (Woo et al., 2007). It has been reported that DDIT3 overexpression can lead to cell cycle arrest and/or apoptosis (Woo et al., 2007). Studies have also shown that DDIT3 can trigger key early events leading to cell apoptosis, which is considered an important target for the development of anti-cancer drugs (Oyadomari and Mori, 2004). Additionally, DDIT3 can participate in cell apoptosis transition and induce Bcl2 down-regulation and DR5 (death receptor 5) activated protein (Farooqi et al., 2015). RASK et al. also indicated that increased DDIT3 was associated with the tumor invasion of CRC (Rask et al., 2000). However, Sun et al. activated the PERK-ATF4-CHOP signaling pathway through TIIA, and then increased the expression of ULBP1 and DR5 through ATF4 and CHOP, leading to enhanced NK cell-mediated killing of NSCLC cells, which seemed to indicate a connection between ULBP1 and CHOP (Sun et al., 2021). In general, the combination of these 3 genes that reflect different levels can improve the prognosis of COAD patients. However, our research still has some unavoidable limitations. Firstly, the study obtained fewer clinical parameters from the TCGA database, and more clinical parameters need to be included to reduce clinical bias. Secondly, we only analyze from the perspective of genes, but due to the limitations of the current experimental conditions, there is no protein-level validation. In the future, more experiments including in vivo and in vitro are needed to explore. Finally, the Guangxi cohort in this study is only a single-center cohort, and multiple centers and larger samples might be needed for further validation.

CONCLUSION
Our study was the first to investigate the diagnostic and prognostic value of the immune-related ULBP1 gene in COAD. ULBP1 gene had a high diagnostic value in COAD. Up-regulated ULBP1 gene of patients with COAD predicted a worse prognosis compared to those patients with downregulated ULBP1 gene. GSEA results showed that ULBP1 was involved in the apoptotic pathway and biological process of T cell mediated cytotoxicity, regulation of natural killer cell activation, and T cell mediated immunity of COAD. The combination survival analysis showed that the combination of high expression of ULBP1, AARS1, and DDIT3 would increase the 2.2-fold death risk of COAD when compared with those of low expression genes. However, these findings need to be further validated.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University. [Ethics No. 2020 (K-Y-E-050)]. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
G-TR wrote the manuscript. G-TR, H-LX, and L-CZ analyzed and interpreted the patient data, G-TR, H-LX, L-CZ, Y-ZG, and H-PS made substantial contributions to the conception, design, and intellectual content of the studies. All authors read and approved the final manuscript.

FUNDING
This work was supported by the National Key Research and Development Program (2017YFC1309200).