HYAL3 as a potential novel marker of BLCA patient prognosis

It has been previously demonstrated that hyaluronan (HA) potentially regulates the initiation and propagation of bladder cancer (BLCA). HYAL3 encodes hyaluronidase and is a potential therapeutic target for BLCA. We aimed to explore the role that HYAL3 plays in BLCA pathogenesis. HYAL3 expression in BLCA specimens was analyzed using The Cancer Genome Atlas (TCGA) database and the Gene Expression Omnibus (GEO) cohort as well as confirmed in cell lines and The Human Protein Atlas. Then, associations between HYAL3 expression and clinicopathological data were analyzed using survival curves and receiver-operating characteristic (ROC) curves. The functions of HYAL3 were further dissected using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and the protein–protein interaction network. Finally, we harnessed the Tumor IMmune Estimation Resource and Gene Expression Profiling Interactive Analysis to obtain correlations between HYAL3 expression, infiltrating immunocytes, and the corresponding immune marker sets. HYAL3 expression varied greatly between many types of cancers. In addition, a higher HYAL3 expression level predicted a poor overall survival (OS) in both TCGA-BLCA and GEO gene chips (P < 0.05). HYAL3 also exhibited an acceptable diagnostic ability for the pathological stage of BLCA (area under the receiver-operating characteristic curve = 0.769). Furthermore, HYAL3 acted as an independent prognostic factor in BLCA patients and correlated with the infiltration of various types of immunocytes, including B cells, CD8+ T cells, cytotoxic cells, T follicular helper cells, and T helper (Th) 2 cells. HYAL3 might serve as a potential biomarker for predicting poor OS in BLCA patients and correlated with immunocyte infiltration in BLCA.

Extracellular matrix (ECM) alteration is closely associated with tumor invasion and progression [6]. In addition, dysregulated ECM is associated with the epithelial-to-mesenchymal transition involving stem cells in cancer. The ECM also regulates tissue metabolism and facilitates the development and progression of different cancers, including BLCA [7]. Hyaluronic acid (HA) is a glycosaminoglycan mainly associated with the ECM. While HA does not induce cellular transformation, it supports other important tumor phenotypes, such as proliferation, migration, resistance to proapoptotic stimuli, and epithelial-to-mesenchymal transition [8].
In previous reports, HA was found to be a reliable tumor marker in patients with urothelial carcinoma [9]. Hyaluronidase (HYAL) is an endogenous glycosidase that degrades HA through the restrictive digestion site at the β-1,4-glycosidic bond between D-glucuronic acid and N-acetylglucosamine. This is achieved by breaking the β-1,4-glycosidic bond between 2-acetyl-D-deoxygenation-D-glucose and D-glucuronic acid [10]. HYAL was first discovered by Duran-Reynals [11]. Through genomic sequencing analyses, researchers have identified multiple HYAL family members, including HYAL1, HYAL2, HYAL3, HYAL4, PH20, and HYALP1. HYAL3 is located on human chromosome 3p21.3 [12,13]. The regulated turnover of HA plays a critical role in many biological processes, including cellular proliferation, migration, and differentiation. Although an association between HYAL3 and tumor development has been reported [14], the mechanism through which HYAL3 regulates tumor phenotypes remains unknown.
According to previous studies, tumor-infiltrating immunocytes, including B, T, CD8 + T cells, and others, play an important role in regulating the balance between antitumor immunity and immune escape in BLCA [15][16][17]. To date, none of the known biomarkers can accurately predict the therapeutic response to immune checkpoint inhibitors in patients with BLCA. However, several reports have shown that cisplatin-based combination chemotherapy might increase CD8 + T cell infiltration and programmed death ligand 1 expression while decreasing the number of immune-suppressing cells [18]. Therefore, it is important to explore the potential mechanisms through which tumor-infiltrating immunocytes regulate the therapeutic responses to chemotherapy and immune checkpoint inhibitors in BLCA.
In the current study, we analyzed the relationships between HYAL3 expression levels, clinical features, and overall survival (OS) in patients with BLCA, utilizing The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), and the Human Protein Atlas databases. This approach was followed by using the Metascape website to enable Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses and the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING)-enabled analysis of the HYAL3-associated protein-protein interaction network. We further used the Tumor IMmune Estimation Resource (TIMER) and the Gene Expression Profiling Interactive Analysis (GEPIA) databases to analyze the associations between the HYAL3 expression level, the type of infiltrating immunocytes, and their corresponding gene marker sets.

Data source
All TCGA expression datasets from RNA-seq were downloaded from the TCGA website (https:// portal. gdc. cancer. gov/). Only cancer datasets including paired samples for an accurate identification of differentially expressed genes were incorporated into our study. No blood samples were included into the analysis as normal samples. A total of 18 cancer datasets matched these criteria. We employed the "RUVg" function in the package RUVSeq (v3.8) to correct the batch effect in the RNA-seq datasets [19]. The identification of differentially expressed genes was performed by the "glmTreat" function in the package edgeR (v3.24.0) [20]. Genes with a count-per-million ≥0.1 in the normal samples were defined as expressed genes. The normalized expression value of trimmed mean of M-values was also generated by the "edgeR" package. We used the above package in R version 3.6.3.

GEO database and the human protein atlas
BLCA-related profiles were obtained from the GEO database (http:// www. ncbi. nlm. nih. gov/ geo/). Data that met the following criteria were selected: (I) studies including at least 20 samples and (II) examination of mRNA expression in both cancer tissue and adjacent normal tissue from BLCA patients. Finally, GSE31684 was selected as our validation cohort. Studies without useful data for analysis were excluded. Differentially expressed genes between BLCA and normal tissue samples were ranked by the Robust Multi-Array Average and Linear Models and annotated by converting the different probe IDs to gene IDs [21].

TIMER and GEPIA databases
TIMER (Version 2.0) was established to explore the abundance of immunocyte infiltration in different tumors. We used TIMER to assess the relationship between HYAL3 and six types of infiltrating immunocytes (CD8 + T cells, CD4 + T cells, B cells, dendritic cells (DCs), macrophages, and neutrophils) in BLCA via gene modules. Expression dispersion maps were created between a pair of custom genes for gastric cancer, and the statistical significance of the correlation and estimation by Spearman analysis were determined by the correlation module. The level of gene expression was shown as log2 RSEM (RNA-Seq by Expectation-Maximization).
GEPIA (http:// gepia. cancer-pku. cn/ index. html) is an online database consisting of more than 8000 types of tumors and normal tissues from the TCGA and the Genotype-Tissue Expression (GTEx) databases. We used these data to explore the association between HYAL3 expression and multiple immunologic marker datasets. The Spearman method was used to determine the correlation coefficient, and the median value of the HYAL3 expression was used as a cutoff to distinguish high expression from low expression.

Cell lines
The human BLCA cell lines T24 and 5637 as well as the noncancerous urothelial cell line SV-HUC were purchased from the Shanghai Institute of Cell Biology, Chinese Academy of Sciences. These cells were cultured in RPMI 1640 (Procell Life Science & Technology, China) and Dulbecco's modified Eagle's medium (Procell Life Science & Technology, China) supplemented with 10% fetal bovine serum (Ausbian Corporation, Australia) and 1% penicillin/streptomycin (Beyotime, Shanghai Biyuntian Biology Technology, China). The cells were maintained at 37 °C in a CO 2 incubator. When the cells reached 80% confluence, the cells were trypsinized and passaged at a 1:3 ratio.

qRT-PCR
Total cellular RNAs were extracted using TRIzol reagent (Invitrogen, USA), according to the manufacturer's instruction. The quantity of RNA was calculated based on the absorbance at 260 nm detected by a Nan-oDrop 2000 spectrophotometer. An absorbance ratio (260 nm/280 nm) between 1.8 and 2.0 was considered as good purity RNA and used for further experiments. Samples of RNA (2 μg) were transcribed into cDNAs with the PrimeScript RT reagent Kit with gDNA Eraser (TAKARA Corporation) in a final volume of 50 μL. Specific cDNAs were amplified with SYBR ® Green Master Mix (TAKARA Corporation) utilizing an ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The reaction conditions were as follows: 95 °C for 2 min; followed by 95 °C for 30 s, and 40 cycles of 95 °C for 5 s and 60 °C for 34 s. The results were analyzed by using the 2 − ΔΔCT relative quantitative method, with GAPDH as an internal control.
Gene-specific primers for HYAL3 and the reference GAPDH were designed using the National Center for Biotechnology Information Primer-Blast Tool (https:// www. ncbi. nlm. nih. gov/ tools/ primer-blast/). All reactions were performed in triplicate, and their melting curves were analyzed to confirm their specificity and accuracy. The expression levels of HYAL3 were normalized to the levels of GAPDH. The sense and antisense primer sequences for HYAL3 and GAPDH were as follows: HYAL3, forward 5′-GGC CAA CGT TGT CGG ACC GAT-3′, reverse 5′-CAG CAT GGC AGC GGC CGG TATAG-3′; and GAPDH, forward 5′-CAG GAG GCA TTG CTG ATG AT3′, reverse 5′-GAA GGC TGG GGC TCA TTT -3′.

Univariate and multivariate cox regression analyses
To determine the influence of HYAL3 on the outcome of patients with BLCA, we used univariate Cox regression analysis to examine the association between HYAL3 and OS in the TCGA-BLCA cohort. Then, we further used multivariate Cox regression analysis to determine whether HYAL3 could independently predict the prognosis of patients with BLCA. The confidence interval (CI) was set at 95%, and a P-value < 0.05 was considered to indicate a significant difference in the statistical analyses.

Survival curves and receiver-operating characteristic (ROC) curves
We constructed Kaplan-Meier survival curves according to the expression levels of HYAL3 to investigate whether HYAL3 expression affected the outcomes of patients with BLCA. Meanwhile, to evaluate the predictability of HYAL3 expression for BLCA prognosis, we incorporated the clinical and pathological data from TCGA-BLCA and HYAL3 expression to generate ROC curves.

Linkedomics database
The Linkedomics Database (http:// www. linke domics. org) includes different sequencing data from 32 cancers in the TCGA and online clinical databases. The HYAL3-related genes were analyzed statistically using Pearson's correlation coefficient, and the data were presented in volcano plots, heat maps, or scatter plots. The results showed the genes exhibiting the closest association with HYAL3 in TCGA-BLCA.

Protein-protein interaction network
We also used the STRING website (https:// string-db. org/) to predict the proteins that interacted with HYAL3. We inputted HYAL3 and then set the confidence score at > 0.4 for significance.

GO and KEGG enrichment analyses
We used the Metascape website (https:// metas cape. org/ gp/ index. html) to conduct GO and KEGG enrichment analyses based on the top 200 genes related to HYAL3 in BLCA. GO consists of three domains: molecular function (MF), cellular component (CC), and biological process (BP). Terms with a P-value< 0.01, a minimum count of 3, and an enrichment factor > 1.5 were collected and grouped into clusters based on their membership similarities.

Statistical analysis
Each experiment was repeated at least three times. Data are expressed as the mean ± standard deviation. Differences between and among groups were compared using the independent-samples T test for qualitative variables. All statistical tests were two-sided, and statistical significance was defined as P < 0.05.

Patient characteristics
The flow chart of the methodologies used in this study is presented in Fig. 1. Our study included data from RNA sequencing and information related to the patient prognosis from 414 BLCA samples and 40 normal tissues in the TCGA and the GTEx databases. All patients were divided into two groups according to the HYAL3 expression levels.

Higher HYAL3 expression levels in tumor samples than in normal tissues
To explore the role of HYAL3 in the development of BLCA, we analyzed the mRNA-sequencing expression levels of HYAL3 in various types of cancer (Fig. 2). The results indicated that HYAL3 might serve as an oncogene in the development of cancers including BLCA. We noted that the expression levels of HYAL3 were significantly higher in BLCA than in normal tissues in the TCGA and the GTEx databases (P = 3.5 e− 09) (Fig. 3a). These findings were also validated by the qRT-PCR assays using the BLCA cell lines (Fig. 3b). Furthermore, we validated the expression levels of HYAL3 in BLCA using the Human Protein Atlas. We found that the expression levels of HYAL3 in BLCA tissues were upregulated compared with those in normal tissues (Table 2). However, when Fig. 6 The univariate (a) and multivariate (b) Cox regression analyses of HYAL3 expression and clinical data the associations between HYAL3 expression and clinicopathological parameters in BLCA patients were compared, we found that the HYAL3 mRNA expression levels were higher in those aged > 70 years old and in those who died compared with the other patients (P = 0.02 and 0.04, respectively; Fig. 4a-b). In contrast, the HYAL3 mRNA levels did not differ according to the TNM stage, histological grade, pathological grade, sex, subtype, or lymphovascular invasion (P > 0.05) (Fig. 4c-j).

Higher HYAL3 mRNA expression correlates with a shorter OS in BLCA patients
To explore the influences of HYAL3 on the OS of BLCA patients, we constructed Kaplan-Meier survival curves according to the expression levels of HYAL3 to investigate whether HYAL3 expression affected the outcomes of patients with BLCA. The results showed that BLCA cases with lower HYAL3 mRNA expression levels had a significantly longer OS in the TCGA-BLCA cohort (P < 0.001; Fig. 5a), and similar findings were also observed in the GSE31684 cohort (P = 0.004; Fig. 5b). According to univariate Cox regression analysis, we found that the TNM stage, age, subtype, pathological stage, and HYAL3 expression level were associated with the OS of the BLCA patients. Multivariate Cox regression analyses indicated that the expression level of HYAL3 could be an independent prognostic factor for BLCA patients (Fig. 6).

HYAL3 as a potential biomarker for predicting the pathological stage of BLCA
We constructed Kaplan-Meier survival curves according to the expression levels of HYAL3 and used ROC curves for quantifying predictive efficacy ( Fig. 7a-g). When the HYAL3 expression was used to predict BLCA presence/ absence, the AUC was 0.647 (95% CI: 0.511-0.783). If the cutoff value was set at 1.238, the sensitivity and specificity were 81.6 and 71.4%, respectively (Fig. 7a). When the HYAL3 expression was used to predict the patient's age, the AUC was 0.559 (95% CI: 0.511-0.783). If the cutoff value was set at 2.525, the sensitivity and specificity were 22.8 and 90.6%, respectively (Fig. 7b). When the HYAL3 expression was used to predict the pathological stage, the AUC was 0.769 (95% CI: 0.629-0.909). If the cutoff value was set at 1.599, the sensitivity and specificity were 57.8 and 99.5%, respectively (Fig. 7c). When the HYAL3 expression was used to predict the T stage, the AUC was 0.596 (95% CI: 0.504-0.587). If the cutoff value was set at 1.57, the sensitivity and specificity were 43.6 and 81.8%, respectively (Fig. 7d). When the HYAL3 expression was used to predict the M stage, the AUC was 0.460 (95% CI: 0.240-0.680). If the cutoff value was set at 1.221, the sensitivity and specificity were 19.8 and 99.5%, respectively (Fig. 7e). When the HYAL3 expression was used to predict the N stage, the AUC was 0.521 (95% CI: 0.469-0.573). If the cutoff value was set at 1.336, the sensitivity and specificity were 63.5 and 43.5%, respectively (Fig. 7f ). When the HYAL3 expression was used to predict the OS, the AUC was 0.566 (95% CI: 0.509-0.622). If the cutoff value was set at 1.826, the sensitivity and specificity were 50.8 and 66.2%, respectively (Fig. 7g). The above results suggest that HYAL3 might be a potential biomarker for predicting the pathological stage of BLCA.

HYAL3 as a tumorigenesis facilitator through interacting with proteins related to immune regulation
To further realize the functions of HYAL3, we used the STRING website (https:// string-db. org/) to predict the proteins that interacted with HYAL3. The top 10 interacting proteins and their gene names, annotations, and scores are listed in Fig. 8, including GUSB, RHCG, ARSB,  Figure 8 presents the interactions between HYAL3 and 10 target proteins with the highest prediction score. SPAG, IDUA, and ELN play vital roles in the regulation of the immune response. These findings indicate that HYAL3 might promote tumorigenesis by regulating the immune system.

HYAL3 potentially regulates tumor differentiation and growth
To explore the downstream target of HYAL3, the Linkedomics Database (http:// www. linke domics. org) was used to predict the targets of HYAL3. The HYAL3-related genes were analyzed statistically using Pearson's correlation coefficient and presented in heat maps. The results showed 200 genes exhibiting the closest association with HYAL3 in TCGA-BLCA (Supplementary Table 1). Heat maps of genes negatively (Fig. 9a) or positively (Fig. 9b) associated with HYAL3 expression are shown in Fig. 10. We used the Metascape website to conduct GO and KEGG enrichment analyses based on the top 200 genes related to HYAL3 in BLCA. The GO enrichment analysis results are listed in Fig. 10. The enriched biological processes mainly involved in tumor growth included chromatin organization, GTP and DNA metabolic processes, cell division, cellular response to UV-B, cellular protein catabolic process, positive regulation of the DNA metabolic process, ribonucleoside diphosphate metabolic process, regulation of the cell cycle process, response to oxidative stress, and the establishment or maintenance of epithelial cell apical/basal polarity. The enriched cellular components mainly included the intracellular protein-containing complex, centrosome, condensed chromosome, preribosome, large subunit precursor, replication fork, nuclear speck, and ficolin-1-rich granule. The enriched molecular functions mainly included catalytic activity acting on nucleic acid, hyaluronoglucosaminidase activity, pyrophosphatase activity, phosphotransferase activity, histone binding, DNA-binding transcription activator activity, ubiquitin-conjugating enzyme binding, exoribonuclease activity, and producing 5′-phosphomonoesters. The  Fig. 11 and include RNA degradation, glycosaminoglycan degradation, biosynthesis of cofactors [22], N-glycan biosynthesis, cell cycle, purine metabolism, and biosynthesis of amino acids. These results indicate that HYAL3 is closely associated with processes related to the regulation of tumor cell differentiation and growth.

HYAL3 expression correlates with immune cell infiltration in BLCA
We analyzed the associations between various types of infiltrating immunocytes and HYAL3 expression in BLCA patients. The results are shown in Fig. 12. To further explore the potential function of HYAL3 in various infiltrating immunocytes in BLCA, we analyzed the GEPIA and TIMER databases regarding the associations between HYAL3 and several immunologic marker sets with the corresponding signs of different types of immunocytes ( Table 3). The results indicate that the HYAL3 expression was associated with Th1, Th2, Th9, Th22 functional T cells, M1 macrophages, and neutrophils in BLCA, suggesting that HYAL3 might be associated with immune responses in BLCA.

Discussion
In this study, we explored the prognostic importance of HYAL3 in BLCA, with regard to the biological processes, molecular functions, cellular components, and potential signaling pathways. We found that the HYAL3 expression level could assist in the diagnosis of BLCA and that it further predicted the OS of BLCA patients. In addition, we explored the correlation between HYAL3 expression and the type of infiltrating immunocytes and found that HYAL3 was associated with glycosaminoglycan degradation. According to the protein-protein interaction network analysis results, it was found that the HYAL3-related gene SPAG9 was associated with immune response [23]. In addition, IDUA is a novel glycolysisrelated gene that is associated with the immune microenvironment in renal cell carcinoma [24]. Moreover, ELN has been reported to be an immune-related gene in BLCA [25]. HYAL3 was associated with several GO functions, including androgen receptor signaling, response to oxidative stress, and negative regulation of the target of rapamycin signaling, the latter of which might correlate with immunocyte infiltration [26][27][28]. In the TCGA-BLCA database, we found that HYAL3 was associated with several types of infiltrating immunocytes, including DCs, mast cells, B cells, T follicular helper cells, interdigitating DCs, activated DCs, regulatory T cells, cytotoxic cells, CD8 + T cells, T cells, central memory T cells, and Th cells. Some of these cells have been shown to participate in the pathogenesis of BLCA [29][30][31]. These results shed light on new pathways through which HYAL3 contributes to BLCA carcinogenesis, potentially by regulating immunocyte infiltration.
BLCA is a high-risk malignancy [32], and an early diagnosis and prognostic prediction are important. Existing prediction models are mostly based on mRNAs, miRNAs, or clinical characteristics [33][34][35]. In addition, infiltrating immunocytes play an important role in regulating BLCA pathogenesis [36]. BLCA has a relatively high tumor mutational burden and is responsive to immunotherapeutic approaches such as Bacillus Calmette-Guerin immunotherapy. Therefore, BLCA is often regarded as an immunogenic tumor [37]. Evidence indicates that a high density of tumor-infiltrating CD8 + T cells is a favorable prognostic factor among BLCA patients, whereas programmed death-ligand 1 expression and tumor-associated macrophages are unfavorable features [38][39][40]. Consequently, the exploration of the association between BLCA and tumor-infiltrating immunocytes might be therapeutically beneficial.
HYAL3 belongs to the group of genes located on chromosome 3p21.3 that exhibit an association with tumor suppression; furthermore, the expression of specific transcript variants may parallel the tumor status [13]. Prior research has suggested that HYAL3 promotes tumor growth in colorectal cancer [41]. A well-coordinated HA  turnover participates in biological processes including cellular proliferation, migration, and adhesion. Moreover, HA may assist in escaping from immune surveillance [14]. HYAL has been shown to accumulate in the tumor microenvironment and to correlate with tumor development and invasion [42]. Researchers also have tried to use HYAL to enhance the antitumor efficacy of tumor vaccines in vivo. HYAL has been demonstrated to potentially increase the permeability of tumor tissues by breaking down HA in the tumor ECM, enabling effective immune responses to be mounted for controlling or eliminating malignant cells [43]. Additionally, during the treatment of glioblastoma multiforme, researchers have found that using HYAL-expressing oncolytic adenovirus ICOVIR17 can degrade HA in tumor cells, subsequently modifying the immunologic landscape of the tumor microenvironment [44]. Altogether, these results support our findings that HYAL3 may participate in encoding HYAL and regulating tumor microenvironment-infiltrating immunocytes. Several limitations existed for this study. For example, the online databases contained multiple issues that warrant consideration, including racial restrictions and the numbers of BLCA cases and normal tissues. A longer follow-up period and an increased number of patients in future studies are still needed.

Conclusion
A higher HYAL3 expression level might predict a shorter OS among BLCA patients. HYAL3 was additionally associated with several types of infiltrating immunocytes in BLCA, including Th cells, T cells, CD8 + T cells, cytotoxic cells, B cells, etc. These data indicate that HYAL3 might serve as a biomarker for BLCA diagnosis and treatment in the future. Nevertheless, more research is needed to verify the biological functions of HYAL3 in BLCA.