Comprehensive Analysis of the Prognosis and Correlations With Immune Infiltration of S100 Protein Family Members in Hepatocellular Carcinoma

S100 protein family members (S100s) are commonly dysregulated in various tumors including hepatocellular carcinoma (HCC). However, the diverse expression, mutation, prognosis and associations with immune infiltration of S100s in HCC have yet to be analyzed. Herein we investigated the roles of S100s in HCC from the Oncomine, Gene Expression Profiling Interactive Analysis (GEPIA), Human Protein Atlas, Kaplan-Meier Plotter, cBioPortal and TIMER databases. Compared with para-cancer tissues, the expression levels of S100A4/S100A6/S100A10/S100A11/S100A13/S100A14/S100P were higher in HCC tissues, while the expression levels of S100A8/S100A9/S100A12 were decreased in tumor tissues. The mRNA levels of S100A2/S100A7/S100A7A/S100A8/S100A9/S100A11 were correlated with advanced tumor stage. Besides, higher mRNA expressions of S100A6/S100A10/S100A11/S100A13/S100A14/S100P were shown to have shorter overall survival (OS), while higher expression of S100A12 was associated with favorable OS. Further, the mutation rate of S100s was investigated, and the high mutation rate (53%) was associated with shorter OS. Additionally, the expressions of S100s were found to be significantly associated with various immune infiltrating cells. Hence, our results showed that S100A6/S100A10/S100A11/S10012/S100A13/S100A14/S100P may be regarded as new prognostic or therapeutic markers and S100s inhibitors may be helpful in the combination of immunotherapies.


INTRODUCTION
Hepatocellular carcinoma (HCC) comprises nearly 80% of primary liver cancer and is the sixth most common and fourth deadly malignancy globally (Bray et al., 2018). Surgical resection or liver transplantation are the only potential measures for radical cure, while these measures are only limited to a small subset of patients. In fact, nearly 75% patients were belonging to middle-late stage of HCC when first diagnosed, and the 5 year survival rate was dismal (Allemani et al., 2015;Katsanos et al., 2017). Thus, novel prognosis or therapeutic biomarkers for HCC are imperative. S100 protein family members (S100s) are well-known for their calcium-binding properties, and widely involved in a variety of physiological and pathological processes, including inflammation, cell proliferation, apoptosis and cancer development. Until now, there are at least twenty members that have been reported (S100A1-S100A14, S100A7A, S100A16, S100B, S100P, S100G, S100Z) (Allgower et al., 2020). The expression of S100s are commonly dysregulated in human malignancies, including HCC (Bresnick et al., 2015;Xue et al., 2015). The expression and role of S100s display a distinctive pattern in different cancer types. For instance, S100A2 plays a tumor-suppress role in most cancers including breast cancer, melanoma, gastric cancer etc., but may have an opposite role in lung cancer (Hountis et al., 2014). S100s have been shown to have complex roles in HCC. Apart from S100A5, S100A7, S100A7A, S100A12, S100A13, S100A16, S100B, S100P, S100G, S100Z, the rest family members, such as S100A1, S100A2, S100A3, S100A4, S100A6, S100A8, S100A9, S100A10, S100A11, S100A14, and S100P are shown to be expressed in HCC (Wang et al., 2001;Kittaka et al., 2008;Hua et al., 2011;Luo et al., 2013;Yan et al., 2013;Zhao et al., 2013;De Ponti et al., 2015;Tao et al., 2017;Guo et al., 2018). Furthermore, the protein expressions of S100A1, S100A4, S100A6, S100A14 and S100P detected by immunohistochemistry were reported to be associated with short survival of HCC patients. Nevertheless, some S100 family members, such as S100A5, S100A7, S100A7A, S100A13, S100A16, S100B, S100G, and S100Z have been rarely studied in HCC. Moreover, the role of S100s in affecting the immune microenvironment in HCC has yet to be explored. However, due to small sample size in some studies, these results remained questionable. As far as we know, bioinformatics analysis has yet to be operated to discover the distinct properties of S100s in HCC. With the help of integrated bioinformatics analysis, we are now able to analyze thousands of interested genes through multiple databases at one time (Liu et al., 2020). On the basis of in-depth analysis, we analyzed the expressions, mutations, predictive signaling pathways and associations with immune infiltration of different S100s in patients with HCC and further explored their possible functions, different prognostic roles and the associations with immune infiltrates in HCC.

Ethics Statement
The study was approved by the Ethics Committee of Zhongshan Hospital. All the data were collected from published databases. It was then confirmed that all written informed consent was obtained.

ONCOMINE Database
Oncomine database, an online database 1 was utilized to analyze the mRNA expression of 20 S100s in different types of cancer including HCC (Rhodes et al., 2004). Students' t-test was applied to get a p-value (cutoff p-value: 0.001, cutoff fold change: 1.5).

GEPIA Dataset
GEPIA, an online tool that contains the RNA expression of 9,736 tumors and 8,587 normal samples from the Cancer Genome Atlas and the Genotype Tissue Expression (GTEx) projects 2 (Tang et al., 2017). The dataset is capable of providing tumor/normal differential expression analysis, survival analysis, correlation analysis and so on. HCC patients were divided into high and low expression groups according to their median value of mRNA expression. The cutoffs of p-value and fold change were defined as following: p-value: 0.01, fold change: 1.

Human Protein Atlas
The Human Protein Atlas, an online tool which contains protein expression data (detected by immunohistochemistry) from 17 different human cancer 3 , was applied to study the protein levels of S100s in HCC (Asplund et al., 2012). Here, we compared the protein expression of different S100s between HCC tissues and normal tissues by immunohistochemistry.

The Kaplan-Meier Plotter
The Kaplan-Meier Plotter, a website that contains gene expression data and survival information of various cancers 4 was used to analyze the prognostic value of different S100s in HCC (Gyorffy et al., 2013). The HCC samples were divided into high and low expression groups according to their median value of mRNA expression and then analyzed by Kaplan-Meier survival plot. Statistically significant was considered when a p < 0.05. cBioPortal cBioPortal, an online tool that contains multidimentional cancer genomics data 5 was used to analyze genetic mutations in S100s and their associations with survival (Gao et al., 2013). The HCC dataset which included 360 patients were chosen for mutation analysis. The co-expression and network was also analyzed.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis
Fifty interactors based on the cBioPortal tool were further explored by GO and KEGG. The functions and predicted pathways were analyzed in the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6 (Huang et al., 2009). We then used the R package "ggplot"and "stringr" to perform GO and KEGG pathway enrichment analyses.
FIGURE 1 | The transcription levels of S100 family members (S100s) in different types of cancers (Oncomine database). Cells in red color represents the overexpression of interested gene, while lowexpression shows in blue. The numbers in these cells represents the numbers of datasets which satisfies the threshold: p-value: 0.001, fold change: 1.5. The levels of S100A4/S100A6/S100A10/S100A11/S100A13/S100A14/S100P were higher in tumor than in normal liver samples, while S100A8 were downregulated.

TIMER Database
TIMER, an online tool 7 that contains comprehensive analysis of immune infiltrates in different types cancers was used to analyze the association between S100s and immune infiltrates in HCC. The associations between gene markers and various immune cells were also recorded . Spearman tests were used for analyzing the correlation coefficients. 7 https://cistrome.shinyapps.io/timer/
Frontiers in Genetics | www.frontiersin.org were significantly upregulated in patients with HCC from multiple datasets, while S100A8 were downregulated (p = 0.001, fold change = 1.5). However, the other members did not differ significantly in HCC as assessed in the multiple datasets (Chen et al., 2002;Wurmbach et al., 2007;Mas et al., 2009;Roessler et al., 2010). We next utilized the GEPIA dataset to analyze the expression profile of twenty S100s. The gene expression profile analysis demonstrated that the levels of S100A4/S100A6/S100A10/S100A11/S100A13/S100P were higher in tumor than in normal liver samples. Moreover, the transcription levels of S100A8, S100A9, and S100A12 were decreased (p = 0.01, Figure 2). Associations between S100s mRNA levels and tumor stage were also analyzed. S100A2/S100A7/S100A7A/S100A8/S100A9/S100A11 groups were shown to be statistically significant (Figure 3). After checking the mRNA levels of S100s in HCC, the Human Protein Atlas was applied to explore the protein expression of S100s in HCC. As displayed in Figure 4, low protein expression of S100A7A was found in normal liver tissues, all the other S100s proteins were not detected in normal liver tissues. S100A5 protein expression was not shown by the database. Low protein expressions of S100A4/S100A11/S100B were observed in HCC tissues. Medium protein expressions of S100A1/S100A6/S100A7/S100A7A/S100A10/S100G/S100P were observed in HCC tissues. Protein expressions of the other members of S100s were not detected in HCC tissues. However, due to the small HCC sample size of immunohistochemistry results in the Human Protein Atlas, the results obtained should be verified further.

Prognostic Roles of S100s in HCC
We explore the prognostic roles of distinct S100s in patients of HCC by using the Kaplan-Meier Plotter tools 8 . As shown in Figure 5 and Table 2, we found that S100A2/S100A6/S100A8/S100A9/S100A10/S100A11/S100A13/S 100A14/S100P high mRNA expressions were significantly associated with poor OS (overall survival). Besides, the high mRNA expressions of S100A5/S100A7/S100A7A/S100A12/S100G/S100Z were correlated with better OS. In addition, the other members 8 http://kmplot.com/analysis/ FIGURE 3 | Correlation between S100s mRNA levels and tumor stage in HCC patients (GEPIA database). The mRNA levels of S100A2/S100A7/S100A7A/ S100A8/S100A9/S100A11 were correlated with advanced tumor stage. A p < 0.05 was considered statistically significant.
Frontiers in Genetics | www.frontiersin.org FIGURE 4 | Representative immunohistochemistry images of distinct S100s in HCC tissues and normal liver tissues (Human Protein Atlas). S100A5 protein expression was not shown by the database. Low protein expression of S100A7A was found in normal liver tissues and medium protein expression was observed in HCC tissues. All the other S100s proteins were not detected in normal liver tissues. Low protein expressions of S100A4/S100A11/S100B were observed in HCC tissues. Medium protein expressions of S100A1/S100A6/S100A7/S100A7A/S100A10/S100G/S100P were observed in HCC tissues.
of the S100s were found to have no correlations with the patient's survival. Based on the S100s mRNA expressions profiles mining above, we concluded that higher mRNA expressions of S100A6/S100A10/S100A11/S100A13/S100A14/S100P were found to be significantly associated with shorter OS, while higher expression of S100A12 was associated with favorable OS.

Genetic Mutations in S100s and Their Associations With Survival of HCC Patients
We further explored genetic alteration in S100s and their associations with survival by using the cBioPortal online tool 9 . As shown in Figure 6A, the mutation rate of S100s was really high in HCC patients. In the 360 enrolled samples, genetic mutations were discovered in 192 HCC patients which indicated a 53% mutation rate. S100A16, S100A13, S100A10 were the top three members with genetic mutations (the mutation rates were 24, 21, and 21%, respectively). Furthermore, high rate of genetic alteration in S100s correlated with shorter OS of HCC patients ( Figure 6B, p = 0.0204). Collectively, the above findings suggested that genetic alteration of S100s may have an impact on prognosis of HCC patients. 9 www.cbioportal.org

Predicted Functions and Pathways of the Mutations in S100s and Their 50 Frequently Altered Neighbor Genes in HCC Patients
We also explored the S100s correlations and networks with the online tool 10 . As was shown in Figure 6C, significant and positive correlations were observed between differentS100s. Positive correlations were shown as red color while negative correlations showed blue. The numbers in the cells represent their correlation coefficient. We next explored the network for S100s and their 50 most associated genes. As shown in Figure 6D, the apoptotic-related genes, including TP53, TP63,TP73, BRCA1, NFKB1, MAPK1, MAPK3, and MAPK14 were correlated with S100s alterations. Moreover, functions and pathways of these genes were further explored by GO and KEGG in DAVID 11 . As was shown in Figure 7A, biological processes, such as innate immune response, positive regulation of transcription, positive regulation of NF-kappaB transcription factor activity, inflammatory response and apoptotic process were remarkably regulated by the S100s mutations in HCC. As was FIGURE 5 | The prognostic value of S100s mRNA expression in HCC patients (Kaplan-Meier Plotter). Higher mRNA expressions of S100A6/S100A10/S100A11/ S100A13/S100A14/S100P were shown to have shorter survival time, while higher expression of S100A12 was associated with longer survival time. OS, overall survival; DFS, disease-free survival. A p < 0.05 was considered statistically significant.
shown in Figure 7B, cellular components, including nucleus, cytoplasm, cytosol, extracellular exosome, and nucleoplasm were regulated by the S100s mutations in HCC. In addition, S100s mutations also regulated the molecular functions, such as protein binding, calcium ion binding, protein homodimerization activity, identical protein binding and transcription regulatory region DNA binding ( Figure 7C). By KEGG analysis, pathways including Toll-like receptor signaling pathway and NF-kappa B signaling pathway (Figures 7D, 8) were found to be associated with S100s alterations.
FIGURE 6 | S100s gene expression and mutation analysis in HCC (cBioPortal) (A) S100s gene expression and mutation analysis in HCC. (B) High genetic alterations (53%) in S100s were associated with shorter OS. (C) Correlation between different S100s in HCC (Positive correlations were shown as red color while negative correlations showed blue. The numbers in the cells represented their correlation coefficient). (D) The network for S100s and the 50 most frequently altered neighbor genes.
FIGURE 7 | The functions of S100s and genes significantly associated with S100s alterations. GO functional enrichment analysis predicted three main functions of S100s mutations and their 50 frequently altered neighbor genes, including biological process, cellular components and molecular functions (A-C). KEGG pathway analysis on S100s and their 50 most frequently altered neighbor genes was shown at (D).
in Supplementary Table 1, significant associations between the S100s and markers of immune cells were found. In short, all the results above suggested that S100s played a critical role in immune infiltration of HCC.

DISCUSSION
S100s dysregulation has been seen in a variety of human malignancies (Allgower et al., 2020). Although the role of some members of S100s in the tumorigenesis and prognosis of several cancers has been partially confirmed, comprehensive understanding the role of twenty S100s in HCC has yet to be conducted. Our work aims to explore the expression, mutation, prognosis and associations with immune infiltration of S100s in HCC. We expect that our results will be helpful in improving existing knowledge, advancing treatment designs, and prognosis prediction. In this study, our results showed that the mRNA levels of S100A4/S100A6/S100A10/S100A11/S100A13/S100A14/S100P were significantly upregulated in HCC patients, while the S100A8/S100A9/S100A12 mRNA levels were decreased. Survival analysis indicated that higher mRNA expressions of S100A6/S100A10/S100A11/S100A13/S100A14/S100P were shown to have shorter OS, while higher expression of S100A12 was associated with longer OS. Moreover, high mutation rate (53%) was found in S100s and correlated with poor survival. These findings suggested that S100A6/S100A10/S100A11/S100A12/S100A13/S100A14/S100P may be used as new prognostic or therapeutic biomarkers. S100A4 is a well-known tumor promoting gene and is the most studied one in HCC. Overexpression of S100A4 led to heavier tumors and more metastasis sites, while down-regulation of S100A4 resulted in a reduction of proliferation and invasion of HCC cells (Yan et al., 2013). S100A4 regulates invasion and migration of HCC cells via NF-kappaB-dependent MMP9 signal (Zhang et al., 2013). It was also reported that S100A4 levels were significantly increased in HCC samples compared with controls,moreover, S100A4 protein levels correlated with tumor differentiation, invasion, recurrence and overall survival . Interestingly, in our study, we showed that both mRNA and protein expression of S100A4 were significantly higher in tumor tissues compared with non-tumor tissues. However, the S100A4 mRNA expression had no correlation with survival time in patients with HCC. One possible reason maybe the protein levels of S100A4 were elaborately controlled by post-transcriptional regulation. And the conclusion obtained here should be verified in further studies. S100A6 is another well-studied member of S100s in cancer. S100A6 is upregulated in many kinds of tumors and may FIGURE 8 | Toll-like receptor signaling pathway and NF-kappa B signaling pathway regulated by the S100s alteration in HCC.
take part in regulating cell proliferation, apoptosis, cytoskeleton dynamics and migration (Bresnick et al., 2015). It was reported that S100A6 was overexpressed in HCC tissues compared with non-tumor tissues and the expression of S100A6 correlated with poor differentiation (Hua et al., 2011). S100A6 enhanced the proliferation and migration of HCC cells via the activation of PI3K/AKT pathway (Li et al., 2014). In our study, the results showed that the expression of S100A6 was significantly upregulated in HCC tissues than in para-tumor tissues. Moreover, high S100A6 mRNA expression correlated with poor prognosis. S100A10 is an oncogenic event in various malignancies, including thyroid anaplastic carcinoma, colorectal cancer and ovarian cancer (Ito et al., 2007;Suzuki and Tanigawara, 2014;Lokman et al., 2016). The role of S100A10 in HCC has rarely been explored. DNA microarray analysis of HCC identified S100A10 as one of the key members linked to integrin and Akt/NF-kappaB signaling (Kittaka et al., 2008). Shan et al. (2013) reported that miR-590-5p inhibited proliferation of HCC cells via down-regulating S100A10 expression. The expression of S100A10 in HCC cells was significantly increased than normal hepatic cells . In our study, FIGURE 9 | Association of 20 S100s genes with immune infiltration in HCC (TIMER). The expressions of S100A2/S100A3/S100A4/S100A5/S100A6/S100A8/S100A9/S100A10/S100A11/S100A13/S100A14/S100A16/S100B/S100P/S100Z were found to be significantly associated with the immune infiltrating levels of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells. A p < 0.05 was considered statistically significant.
the results showed that the expression of S100A10 was significantly upregulated in HCC tissues than in para-tumor tissues. Moreover, high S100A10 mRNA expression correlated with poor prognosis. S100A11 has dual roles in cell growth regulation, intracellular S100A11 inhibits cell growth and extracellular S100A11 promotes cell proliferation (Sakaguchi and Huh, 2011). S100A11 is reported to be associated with advanced tumor behavior and poor survival in bladder cancer (Yao et al., 2007). S100A11 activates Wnt/beta-catenin Signaling and promotes proliferation and migration of cervical cancer cells (Meng et al., 2019). In HCC, S100A11 was found to be a metastasis-associated protein by proteomic analysis (Song et al., 2006). Overexpression of S100A11 promotes migration of HCC cells (Luo et al., 2013). A recent study showed that S100A11 overexpression promoted inflammation/fibrosis and was associated with tumor grade and poor survival time of HCC patients (Sobolewski et al., 2020). In our study, the results showed that the expression of S100A11 was significantly upregulated in HCC tissues than in para-tumor tissues. Moreover, high S100A11 mRNA expression correlated with advanced tumor stage and poor prognosis. S100A12 is commonly expressed on neutrophils monocytes and macrophages, the plasma levels of S100A12 are systemic biomarkers for many inflammatory diseases (Meijer et al., 2012). S100A12 is dysregulated in many types of cancers, upregulated in colorectal cancer and papillary thyroid cancer but reduced in gastric carcinoma and oropharyngeal squamous cells (Thierolf et al., 2008;Funk et al., 2015;Li et al., 2016;Wang et al., 2020). Studies on S100A12 in HCC are limited. Proteomics analysis indicated that S100A12 might be useful marker for predicting the early recurrence/metastasis of HCC after resection (Huang et al., 2014). Another study revealed that S100A12 was expressed not on HCC cells but intratumoral stroma cells and S100A12 overexpression predicted poor prognosis (Cai et al., 2018). In this study, the results showed that the mRNA expression of S100A12 in HCC tissues was significantly lower than that in normal tissues. Moreover, high S100A12 expression correlated with favorable OS in patients with HCC.
There are limited studies on S100A13 in cancer. S100A13 was first identified as a powerfully angiogenic biomarker in some kinds of tumor such as melanoma and gliomas (Landriscina et al., 2006;Massi et al., 2010). A recent study showed that S100A13 overexpression promoted migration of lung cancer cells (Pierce et al., 2008). Consistent with these results, the results showed that the expression of S100A13 was significantly upregulated in HCC tissues than in para-tumor tissues. Moreover, high S100A13 mRNA expression correlated with poor prognosis. S100A14, a novel member of the S100s, has been characterized by differentially expressed and functioned in various cancer types. Overexpression of S100A14 has been found in some kinds of cancers, such as breast cancer, lung cancer and bladder cancer. But downregulated in some other kinds of cancers, for example colorectal carcinoma and renal cancer (Basnet et al., 2019). S100A14 activates ERK and NF-κB signaling and promotes ESCC cell proliferation. However, in oral squamous cell carcinoma, S100A14 positively affects the function of p53 and inhibits cell proliferation (Jin et al., 2011;Sapkota et al., 2012). Few studies investigated the role of S100A14 in HCC. Zhao et al. (2013) reported that knock-down of S100A14 decreased proliferation and invasion of HCC cells. Similar to these results, the results showed that the expression of S100A14 was significantly upregulated in HCC tissues than in para-tumor tissues. Moreover, high S100A14 mRNA expression correlated with poor prognosis. S100P is an oncogenic event in many types of tumors, including pancreatic adenocarcinoma, breast cancer and ovarian cancer (Wang et al., 2006;Ezzat et al., 2016;Bai et al., 2018). In ovarian cancer, S100P was shown to be positively correlated with CA125 and poor prognosis . Studies investigating the role of S100P in HCC were limited. S100P was found to be overexpressed in HCC tissues and correlated with advanced tumor behavior and poor prognosis (Yuan et al., 2013). Another study showed that S100P was useful predictor for microvascular invasion and portal vein tumor thrombus (Qi et al., 2021). Similar to these results, our findings showed that the expression of S100P was significantly upregulated in HCC tissues than in para-tumor tissues. Moreover, high S100P mRNA expression correlated with poor prognosis.
HCC is now generally considered as an immunogenic tumor (Pardee and Butterfield, 2012). Various immune cells have been proved to affect the progression and prognosis of HCC. For instance, tumor-associated macrophages (TAMs) and tumor infiltrations of neutrophils are associated with poor prognosis in patients with HCC, while CD8 + cytotoxic T lymphocytes (CTLs) and dendritic cells (DCs) were predictive for better prognosis (Chen et al., 2020). In this study, our results indicated that S100s played a critical role in immune infiltration of HCC. Especially, both S100A4 and S100A6 had strong positive correlations with macrophages. Therefore, the combination these S100s inhibitors with immunotherapies may be a promising strategy in HCC therapy.
Nevertheless, our study had some limitations. First, all the data included for analysis were retrieved from online databases, further researches are still required to verify these findings. Second, all survival analyses were conducted by the Kaplan-Meier database, the independent prognostic values of S100s by multivariate analysis were not acquired because of data limitation. Third, exploring the detailed mechanism of each S100s in HCC was not conducted in this research, future researches are still needed.
In conclusion, in this study, we systemically explored the expression, mutations, predictive signaling pathways, prognosis and associations with immune infiltration of different S100s in patients with HCC. Our research gives new insights regarding the roles of S100s to HCC and may be helpful for the further combination of immunotherapies with S100s inhibitors in HCC.

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
Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.