Identification and Validation of a Stromal EMT Related LncRNA Signature as a Potential Marker to Predict Bladder Cancer Outcome

Bladder cancer (BLCA) has become one of the most common malignant tumors in the genitourinary system. BLCA is one of the tumors considered suitable for immunotherapy because of the large proportion of immune cells in TME. Epithelial to mesenchymal transition (EMT) is closely related to tumor immunity through its crosstalk with immune cells. A recent study validated that EMT-related genes were mainly expressed by stromal cells and could influence immunotherapy responsiveness. Stromal EMT-related gene signature was also demonstrated to affect the prognosis of multiple tumors, including BLCA. To further explore the prognostic roles of stromal components, we performed a comprehensive analysis of LncRNAs closely associated with stromal EMT-related genes in the TCGA BLCA cohort. We identified a signature including five stromal EMT gene-related LncRNAs that showed significant prognostic value for BLCA patients. By the CIBERSORT and MCP-COUNTER algorithm, we found the signature was markedly correlated with infiltrated immune cells and stromal components of the tumor microenvironment, which may further influence patient’s responsiveness to immune checkpoint blockade therapy. Through immunohistochemical analysis, we confirmed the correlation of the signature with macrophages M2 and CAFs. Meanwhile, key genes related to these LncRNAs, including VIM, MMP2, were also differentially expressed in the stromal components concerning the signature. Our research confirmed the prognostic and immune-associated role of stromal EMT-related LncRNAs. Meantime, we further confirmed that EMT-related genes were mainly expressed in stromal components. Targeting these LncRNAs as well as their related stromal EMT genes may provide potential therapeutic targets for BLCA immunotherapy and precision medicine.


INTRODUCTION
Bladder cancer (BLCA) has become one of the most common malignant tumors in the genitourinary system (1). Transitional cell carcinoma is the most common pathological type of BLCA. According to the invasion depth, BLCA can be divided into nonmuscle invasive(NMIBC) and muscle invasive bladder cancer (MIBC).NMIBC can be treated by transurethral resection of bladder tumor (TURBT) and has a favorable prognosis but still faces recurrence risk. MIBC is poorly treated clinically and often with radical surgery supplemented by postoperative chemotherapy and radiotherapy (2). The 5-year survival rate of MIBC patients is only 50% (3). Due to the rapid development of immunotherapy, immune checkpoint inhibitors for advanced BLCA have recently received significant attention. However, the responsiveness to existing immune checkpoint inhibitors is limited among BLCA patients, partly due to a complex heterogeneous tumor microenvironment (TME). Further study of TME is of great importance for BLCA immunotherapy. TME is a complex and integrated system mainly composed of stromal cells and infiltrated immune cells. Emerging evidence suggests that the stromal components can shape antitumor immunity and affect immunotherapy responsiveness, thus promoting tumors' malignant development (4). The TME is particularly important in BLCA because of the overwhelming evidence that BLCA represents a growing number of solid tumors characterized by a significant number of stromal and immune cells in the TME (5,6). The heterogeneity of TME in BLCA is closely related to patients' different response rates to immunotherapy, in which expression of stromal epithelial to mesenchymal transition (EMT) related genes plays a vital role (7).
EMT has been defined as a dynamical process with intermediate states (8). A complicated regulatory network is engaged in EMT's dynamic procedure at different levels, which further remodels the tumor extracellular matrix and promotes tumor metastasis. Although EMT is a biological process unique to tumor cells and has been studied extensively in vitro and in model organisms, evidence for this phenomenon in human tumors has been limited (9). In contrast, accumulating evidence showed that the up-regulation of EMT-related genes in bulk tumors was driven by expression changes in fibroblasts rather than in epithelial tumor cells (10). For example, Li et al. demonstrated that EMT-related genes were up-regulated only in a subpopulation of cancer-associated fibroblasts (CAFs) through single-cell RNA sequencing (11). A recent study also indicated that stromal EMT-related gene expression might alter T-cell infiltration level, which might eventually impact the responsiveness to immune therapy and patient's survival (7). These findings illuminated intertwined and complicated crosstalks between the stromal components and the EMT process in the development of cancer immune evasion, making the stromal components promising therapeutic targets for cancer immunotherapy.
LncRNAs are novel, potential therapeutic targets and biomarkers for cancer treatments (12). Moreover, researchers have demonstrated that LncRNAs obtain more specificity on indicating actual tumor condition than other types of markers (13). Previous studies have already discussed the prognostic value of stromal EMT-related genes in BLCA (7,14), but the roles of stromal EMT-associated LncRNAs (sEMTLncs) have rarely been reported. Therefore, we carried out this present study to explore the function of sEMTLncs and seek their potential application for predicting BLCA outcome. We constructed a stromal EMT-related LncRNA prognostic signature (sEMTLS) in the present study, which showed good predictive accuracy for BLCA outcome. By adopting the CIBERSORT and MCP-COUNTER algorithm, we further found the sEMTLS was significantly associated with the levels of CAFs and tumor infiltrated immune cells (TIICs). Through immunohistochemical validation of the expressions of ACTA2, CD206, and these sEMTLnc-targeted EMT genes, we found that with the different expression levels of these sEMTLncs, these genes were also differentially expressed in the stromal components. These results shed light on the dual regulatory effect of sEMTLncs on both stromal and immune components in BLCA and laterally confirmed the EMT-related genes' expression was indeed in tumor stromal but not in epithelial tumor cells. Moreover, the ImmuneCell AI database (15) predicted that the expression of these sEMTLncs might also affect BLCA patients' responsiveness to immune checkpoint blockade (ICB) treatment. Targeting these sEMTLncs as well as their related stromal EMT genes may provide potential therapeutic targets for BLCA immunotherapy and precision medicine.

METHODS AND MATERIALS
Raw Data Acquisition 430 samples of transcriptome profiling data, including 19 normal samples, 411 tumor samples, and the corresponding clinical data of 405 bladder transitional cell papilloma and carcinoma patients were downloaded from the TCGA database (https://portal.gdc. cancer.gov/). Finally, 403 BLCA patients were selected and randomly arranged into training (n = 203) and testing groups (n = 200) for further study. Patients' baseline information was listed in Table 1.

Stromal Scoring and Differential
Expressed Genes (DEGs) Screening R language version 4.0.2 loaded with ESTIMATE package was used to calculate the scores of the immune and stromal component in TME for each sample of BLCA patients. After estimation, stromal score, immune score, and estimate score were obtained, representing the abundance of stromal, immune components, and the total. Then, samples were categorized into high and lowstromal score groups based on all samples' median stromal score. Package limma was applied to perform differential analysis of the gene expression; DEGs were screened by comparing the gene expressions in samples between high and the low stromal score groups. DEGs with fold change more than 1 after transformation of log2 and FDR <0.01 were considered significant.
sEMTLnc Signatures sEMTLncs affecting the survival of BLCA were selected by univariate COX analysis using R software survival packages (P < 0.01). Lasso and multivariate cox regression were further used to construct the sEMTLS. Hazard ratio (HR) was used to classify sEMTLnc into the protective (HR < 1) and deleterious (HR > 1) portion ( Table 2). Patients were classified into high-and low-risk groups based on the medium riskscore of the signature. The riskscore was calculated by the formula as followed:

Survival Analysis
R language v4.0.2 with package survival and survminer were used for survival analysis. Kaplan-Meier survival analysis was used for analyzing the survival difference between different groups. P-value of the log-rank test less than 0.05 were considered significant.

Time and Multi-ROC Curves
ROC curves analyzed the predictive accuracy of the signature on BLCA overall survival at 1, 3 and 5 years. Package timeROC of R language v4.0.2 was used for plotting timeROC curves. Independent risk analysis of different clinical characteristics, including age, gender, stage, T classification, and risk score on predicting 1-year overall survival, was conducted by Package survivalROC ( Table 3).

Principal Component Analysis (PCA) and Nomogram Construction
PCA was used to cluster the samples based on the expression of sEMTLnc. A 3D scatterplot visualized patients' distribution. Nomogram was constructed by including the expression level of the sEMTLnc in the signature. A calibration plot was used to explore the calibration and discrimination of the nomogram.

GO, KEGG, and GSEA Enrichment Analysis
DEGs between the high-and low-risk groups of sEMTLS were used for GO and KEGG enrichment analysis. The  Hallmark_Epithelial_Mesenchymal_Transition gene set was used for GESA between the low-risk and high-risk group, which was performed using the GSEA software (version 4.1.0) obtained from the Broad Institute. NOM p <0.05 and False Discovery Rate (FDR) q <0.05 were considered to be significant.

Calculation of TIIC Levels and CAFs Abundance
The CIBERSORT algorithm was used for calculating the levels of 22 different types of TIICs in all tumor samples; Samples with a p-value less than 0.05 were selected for the following analysis. The MCP-COUNTER algorithm (16,17), provided by TIMER 2.0 (http://timer.cistrome.org/) (18), was applied for calculating the abundance of TME components, including endothelial cells and CAFs.

ICB Treatment Reactiveness Predicting
Immune cell abundance identifier (Immune cell AI, http:// bioinfo.life.hust.edu.cn/ImmuCellAI) was applied to estimate the difference of immune cell infiltration in low-and high-risk groups. Patients' response to ICB therapy was then predicted based on the estimation of immune cell infiltration levels.

Real-Time Quantitative PCR
According to the manufacturer's instructions, we used triazole (Invitrogen) to extract total RNA from all recruited BLCA samples. cDNA Synthesis Kit (Osaka, Japan of TaKaRa) combining with RNA (1 mg) was utilized to reverse-transcribed cDNA. The quantitative polymerase chain reaction (qPCR), using the SYBR-Green method (TaKaRa), was performed on an ABI 7500 real-time PCR system (Applied Biosystems). The relative expression level of each lncRNA was calculated by the 2 −DDCt method after normalizing to b-actin level. The forward and reverse primer sequences are shown in Table 4.

Immunohistochemistry Analysis
The gene expression in tumor tissues was detected using the BenchMark GX automatic multifunctional immunohistochemical staining system (Roche, Switzerland) with OptiView DAB Detection Kit (Ventana, USA) according to the manufacturer's instructions. The staining's straightforward procedures were listed as follows: deparaffinization and epitope retrieval in cell conditioner for 90 min. Short (8 min), mild (30 min), and standard (60 min) cell conditioning was performed after epitope retrieval. Primary antibodies were then incubated with the section for 32 min followed by biotinylated anti-IgG antibody and streptavidin-  biotinylated-complex horseradish peroxidase. Hematoxylin was used for counterstaining and Bluing Reagent for post counterstaining. Details of the analyzed genes were listed in Table 5.

Statistics Analysis
Correlation of sEMTLS with infiltrated immune cells was analyzed with Pearson correlation test; p <0.01 was considered significant. Kaplan-Meier curve with log-rank test was used to evaluate the OS between different groups. The Wilcoxon test examined the differences for variables of two groups. Kruskal test estimated statistical significance for variables of more than two groups. Univariate and multivariate Cox regression analyses were displayed to verify the independent prognostic factors for BLCA. Fisher exact test was used to calculate the difference of ICB response between high-and low-risk groups. Two-sided Pvalue <0.05 was considered significant. R language v4.0.2 was used for all statistical analyses.

sEMTLnc Were Identified for Signature Construction
After we scored the stromal components by the ESTIMATE algorithm, DEGs between low and high stromal score groups were screened. 2,693 DEGs were identified. The top 50 of both up-and down-regulated genes with the most significant fold changes were represented by the heatmap ( Figure 1A). The distribution of all DEGs on the two dimensions of -log10(FDR) and logFC was depicted in the volcano plot ( Figure 1B). 421 EMTLncs were enrolled by co-expression with EMT-related genes. 82 sEMTLncs were selected after the intersection of the EMTLnc with the DEGs ( Figure 1C). After the intersection with ImmLnc, 72 out of 82 sEMTLncs were confirmed to be highly coexpressed with genes related to the immune process, which further demonstrated the close relationships between EMT and the immune process. These 72 sEMTLncs were chosen for the construction of sEMTLnc signature.  Table 2). A forest plot illustrated the corresponding HRs and 95% CIs for each sEMTLnc ( Figure  2A). Patients were classified into a high-risk group and a low-risk group based on the training group's median risk score. Patients' overall survival (OS) in the high-risk group was significantly shorter than that in the low-risk group (p < 0.001). Time ROC curve showed good predictive accuracy with AUC of 0.777, 0.776, and 0.799 for predicting 1, 3 and 5 years' OS. Subsequently, we validated the sEMTLS in testing and the entire groups. Statistically, significant OS differences were observed between the high-and low-risk groups in the testing (p < 0.001) and the entire group (p < 0.001). AUC of time ROC curves in the testing and the entire group were 0.667, 0.666, 0.663, and 0.709, 0.719, 0.740 for predicting 1, 3, and 5 years' OS, respectively ( Figure 2B).

The Risk Score of sEMTLS Could Essentially Predict the Clinical Status of BLCA Patients
Based on our signature, the mortality risk of patients in each group climbed with the risk score increased. The expression levels of AL583785.1 and LINC01711 were elevated, while TMEM51-AS1, AC073534.1, and LINC02446 expressed decreasingly as risk score increased ( Figure 3A). Combined with other clinical and demographic characteristics of BLCA patients, the risk score was identified by multivariate cox regression analysis to be an independent prognostic factor for BLCA patients, with a hazard ratio of 1  Figure 3B). The ROC curve validated the prognostic accuracy of the risk score. The AUC of the risk score was higher than any other clinical and demographic characteristics in each group, which further suggested the risk score could be an independent prognostic factor ( Table 3). To further validate the prognostic value and explore the broad applicability of sEMTLS, we analyzed the relationships between sEMTLS with different clinical features, including patient age, tumor grade, staging, and TNM classification in the entire group. The risk score was found significantly related to all of these clinical features confirming the significant association of sEMTLS with the progression of BLCA (Figures 4A-F). Last, we constructed a nomogram of sEMTLS to predict patient survival ( Figure 4G). The calibration curve indicated that

PCA and Functional Analysis Between High-and Low-Risk Groups of sEMTLS
Based on the expressions of sEMTLnc recruited in the sEMTLS, we employed the principal component analysis (PCA) to provide an overview of different distribution patterns between the lowrisk group and the high-risk group. Results indicated a relatively scattered distribution of patients in testing and the entire groups, but with a small overlap in the training group ( Figures 5A-C). GO and KEGG analysis of DEGs between high and low-risk groups of the entire group were enriched in terms related to extracellular matrix remodeling, including extracellular matrix organization, extracellular matrix structural constituent, and ECM-receptor interaction ( Figure 5D). GSEA further proved the functional annotation, with the more EMT-related activity in the high-risk group ( Figure 5E).

Close Relationships Were Found Between sEMTLS and TIICs
Utilizing the CIBERSORT algorithm, we obtained an estimation of the abundances of 22 TIICs. The infiltration proportion of the immune cells in each sample was shown in the barplot ( Figure  6A). Among all the TIICs, Macrophages M0 and M2 were positively correlated to the risk score while T cell CD8+ and T cell CD4 memory activated exhibited negative correlations ( Figure 6B). Further, elevated levels of Macrophages M0, M2, and decreased levels of T cell CD8+ and T cell CD4 memory activated were found in the high-risk group when compared with the low-risk group ( Figure 6C). Combined with the survival time and survival state, all these four TIICs showed significant relations with the OS of BLCA patients, with Macrophages M0, M2 being detrimental factors and T cell CD8+, T cell CD4 memory activated being protective factors ( Figure 6D). The above results highlighted the association between sEMTLS and TIICs, indicating the immune-modulating role of the sEMTLnc.

Correlation of sEMTLS With CAFs and Its Predicting Value to the Responsiveness of ICB Treatment
MCP-COUNTER algorithms calculated CAFs abundance in patients of the TCGA BLCA cohort. The relative abundance of CAFs was represented in the heatmap ( Figure 7A). After comparing the TME components calculated by MCP-COUNTER between high-and low-risk groups, T cells CD8+ were further confirmed to be reduced in the high-risk group, while CAFs abundance was significantly higher in the high-risk group than in the low-risk group ( Figure 7B). The risk score was further validated positively correlated with CAFs abundance and the stromal score ( Figure 7C). Using the immune cell AI database, We found significantly lower ICB treatment responsiveness in high-risk groups than the low-risk group (p=0.040). The risk scores between responders and non-responders also differed with relative higher risk scores in non-responders (p = 0.014) ( Figure 7D).

Validation of the Association Between sEMTLS and TME Components in a Clinical Cohort
A clinical BLCA cohort of 16 patients with different stages was established to validate the correlation between sEMTLS, sEMTLnc targeted key genes, TIICs and CAFs. In our study, a total of five molecules were used for further research, of which CD206 and ACTA2 were used to represent the relative expression of macrophage M2 and CAFs, while three key lncRNA-targeted genes, VIM, CALU and MMP2, were used to detect expression in patients at different risks. The above five molecules in the TCGA cohort were significantly differentially expressed between the high-and low-risk groups ( Figure 8A). Relative expressions of the five sEMTLnc were analyzed by realtime quantitative PCR. The risk score of each patient was calculated according to the formula. A close relation of sEMTLnc with the clinical stage was found and shown in the barplot ( Figure 8B). Expression of CD206 was found in areas where ACTA2 was expressed in high risk patients ( Figure 8C). We further compared the expression of CD206, ACTA2, VIM, CALU and MMP2 between patients with low and high risk scores ( Figures 8D-H). The results confirmed an elevated expression of all the five genes in high risk group patients. Furthermore, the expressions of ACTA2, VIM, and MMP2, which are markers of CAFs, were mainly found in the stromal area. These results demonstrated the relationship between CAFs and macrophages M2 and highlighted the association of sEMTLS with CAFs and macrophages M2 in BLCA patients. Our immunohistochemical results also revealed that in addition to CAFs, vascular smooth muscle cells also expressed EMT-related genes, which confirmed and complemented the previous findings of Li et al.

DISCUSSION
BLCA is one of the tumors considered suitable for immunotherapy because of the large proportion of immune cells in TME. Recent advances in treatment indeed demonstrated immune checkpoint blockers on T cells result in significantly improved survival in BLCA (19). However, challenges remain since there are still many BLCA patients who showed low responsiveness to ICB therapy. Over the past decade, EMT has been considered as a pivotal regulator in metastatic progression (20) and therapy resistance (21), including chemo-(22), radio-(23), and targeted therapy resistance. Under the current understanding, the EMT process is profound for regulating immune cells in TME, including CD8+ T cells and macrophages M2 (24). Evidence is now accumulating that such crosstalks might be a critical mechanism in promoting cancer immune escape (25).
Recent studies have identified EMT as a dynamic process with an intermediate status, making the EMT process a promising target for therapeutic intervention (26), which may further promote immunotherapy efficacy in low responsive cancer patients (27). The TME has been characterized as inflammatory and immunosuppressive (28,29). It owned a variety types of immune and stromal cells. Among the immune cells, macrophages are highly plastic and crucial to the EMT process (30). Previous studies have confirmed that M2 macrophages secrete a series of cytokines that promote EMT and cancer progression via multiple signaling pathways, including ZEB1 (31), SNAL1 (32), VIM (33), and TWIST1 (34). Stromal cells dominated by CAFs also participate in the EMT process. One previous study demonstrated that CAFs could increase the expressions of EMT-related genes and differentiated the recruited monocytes into M2 macrophages, further exerting their immunosuppressive roles targeting CD8+ T cells (35). These results are consistent with our finding in this study that macrophage M2 and CAFs are co-expressed in the stromal compartment. Given all the above evidence, we clearly noticed the association between TME components and the EMT process. Although EMT is a biological process unique to tumor cells has been studied extensively in vitro and in model organisms, evidence for this phenomenon in human tumors has been limited (11). A recent study indicated that stromal components of BLCA comprise a key source of EMT-related gene expression. Expression of these stromal genes correlated with T-cell infiltration and impacted response to immune checkpoint blockade, further influencing BLCA patients' survival (7). In the present study, we found that the expression of sEMTLnc may also affect the infiltration of T-cells and macrophages M2. Immunohistochemical results of these lncRNA-targeted EMT genes also validated that EMT genes were mainly expressed in stromal components. Whereas, questions may raise that since EMT is a unique process of the tumor cells, why EMT related genes were mainly expressed in stromal components but not in tumor cells themselves. Available evidence suggests that the EMT process may not only act in the epithelial tumor cells but also involve intertwined and complex crosstalks between stromal components and tumor cells. Based on the fact that EMT-related genes are only expressed in a subpopulation of CAFs, CAFs have been suggested to even arise from tumor cells undergoing EMT (7). It is difficult to explain these mechanisms by bulk RNA sequencing. Further research is still needed to validate the exact mechanisms between stromal EMT-related gene expression and tumor cells' EMT process. Perhaps applying single-cell sequencing to explore these EMT-associated genes' cellular origin could be useful in explaining this intricate process.
On the other hand, LncRNAs are novel, potential therapeutic targets and biomarkers for cancer treatments (12). Moreover, researchers have demonstrated that LncRNAs obtain more specificity on indicating actual tumor condition than other types of markers (13). Previous research indicated that stromal EMTrelated genes could provide a predictive role in tumor patients' prognosis and affect response to ICB therapy in BLCA patients. Could sEMTLnc play a similar role? From the present results, we confirmed the prognostic value of sEMTLnc in BLCA patients. The sEMTLS could be an independent risk factor to patients' OS. Close correlations between sEMTLS, macrophage M2 and CAFs were also found, suggesting that the combined expression of sEMTLnc was associated with the abundance of macrophages and CAFs. Also, there was a negative correlation between sEMTLS and CD8+ T cells and a significant difference in ICB treatment responsiveness between high and low-risk groups. These results suggested that sEMTLnc, like stromal EMT-related genes, had a significant influence on the immunotherapeutic response and may ultimately affect the prognosis of BLCA. Using a clinical validation cohort, we further confirmed the relationship between CAFs, macrophage M2 and the sEMTLnc.
The key genes related to the sEMTlnc, including VIM, MMP2, and CALU, also expressed differently concerning the risk score. VIM and MMP2 are also markers of CAFs, which promote cancer progression and metastasis (36,37). Simultaneously, CALU was demonstrated to express significantly higher levels in the metastatic cancer tissues (38). Through our immunohistochemical results, we could see that markers such as ACTA2, VIM, and MMP2 were mostly differentially expressed in the stromal area, including CAFs and vascular smooth muscle cells, but not in the tumor cells between the different risk score groups. These results are consistent with the previous literature (7,11), while the detection of expression of key genes targeted by these sEMTlnc in vascular smooth muscle cells also suggested that these EMTlnc may be involved in the formation of tumor neovascularization (39). Taken together, we suggest that these sEMTlnc play a key role in TME remodeling and regulation of TIICs, including T cells and macrophages. The expression of sEMTlnc may ultimately influence immune responsiveness and overall survival of BLCA patients. Further studies on sEMTLnc may provide potential therapeutic targets for BLCA immunotherapy and precision medicine. However, limitations still existed since our results were mainly based on the bioinformatics analysis. Although highly correlated with EMT-related genes, critical questions still existed: Do sEMTLnc expression indeed reflect the biological process of EMT, How do sEMTLnc expression alter T-cells and macrophages M2 infiltration, why sEMTLnc can impact outcomes and affect the responsiveness to ICB treatment in BLCA patients. Simultaneously, It is not sufficient to study sEMTLnc-associated TIICs only by bulk RNA sequencing data based on bioinformatics algorithms, which may lead to different results between different algorithms. The exact association between sEMTLnc, CAFs, macrophages, and T cells still needs to be verified by a series of in vitro and in vivo experiments, including single-cell RNA sequencing. Besides, due to the limited number of cases in the clinical cohort we included, we still need a larger external clinical cohort to validate the relative expressions of the sEMTLnc and the predictive value of the signature.

CONCLUSION
In the present study, we constructed a prognostic signature containing five sEMTLnc, which predicted BLCA patients' prognosis. Further study on the signature confirmed its significant correlation with the abundance of CAFs, macrophages M2 and CD8+ T cells. Similar to the literature reporting that stromal EMT gene expression affects BLCA prognosis and immunotherapy responsiveness, the combined expression of stromal EMT-related lncRNAs may also affect BLCA outcome and immunotherapy responsiveness. Through immunohistochemical analysis, we laterally verified that EMTrelated genes are mainly expressed in the stromal components, including CAFs and vascular smooth muscle cells. The crosstalk between tumor stroma and tumor cell's EMT process is intricate and requires in-depth study. Further study of stromal EMT-related lncRNAs and their targeted genes will help provide possible new targets for BLCA precision therapy and immunotherapy.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Medical Ethics Committee of Shanghai General Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
YD and BW had an equal contribution to this manuscript. HL designed the whole study. YD and BW participated in the bioinformatics and statistical analysis, XJ and JC did the immunohistochemistry analysis. YD performed real-time PCR analysis. YW and JY made the manuscript and figure editing. XW and HL revised the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was funded by the National Natural Science Foundation of China (Grant number: 81972371) and Basic Research on medical and health Application of Suzhou Municipal Science and Technology Bureau (Grant number: SYSD2020076).