Pan-Cancer Analysis and Experimental Validation of SOX4 as a Potential Diagnosis, Prognosis, and Immunotherapy Biomarker

Simple Summary Data show that the global cancer burden continues to grow. It is estimated that there were 19.3 million new cancer cases and nearly 10 million cancer deaths worldwide in 2020. Therefore, a comprehensive analysis is of great significance for the accurate diagnosis and effective treatment of tumors. This study aims to visualize the expression, survival, mutation, methylation, ceRNA network, and immune and prognostic models, as well as to explore the potential role of SOX4 in different tumor types and specific LIHCs using different online tools. In addition, we also confirm that SOX4 may be related to lenvatinib resistance in HCC patients through in vitro experiments. Knockdown of SOX4 can significantly inhibit the proliferation of HCC cells. Finally, this study highlights the key roles of SOX4 in the diagnosis and prognosis of tumors, especially in LIHC, and its potential as a promising therapeutic target for tumors. Abstract Introduction: SOX4 plays an important role in tumorigenesis and cancer progression. The role of SOX4 in pan-cancer and its underlying molecular mechanism in liver hepatocellular carcinoma (LIHC) are not fully understood. In this study, a comprehensive analysis and experimental validation were performed to explore the function of SOX4 across tumor types. Methods: Raw data in regard to SOX4 expression in malignant tumors were downloaded from the TCGA and GTEx databases. The expression levels, prognostic values, genetic mutation, and DNA promoter methylation of SOX4 across tumor types were explored via systematic bioinformatics analysis. The ceRNA regulatory network, immune characteristics, and prognostic models were analyzed in LIHC. Finally, we conducted in vitro experiments including Western blotting, cell proliferative assay, trypan blue staining, and fluorescence microscopy to further explore the function of SOX4 in LIHC. Results: SOX4 expression was significantly upregulated in 24 tumor types. SOX4 expression level was strongly associated with unfavorable prognoses, genetic mutations, and DNA methylation levels across different tumor types. Especially in LIHC, LINC00152/hsa-miR-139-3p/SOX4 was identified as a crucial ceRNA network. Moreover, this study also provides insight into the roles of SOX4 expression in immune cell infiltration, macrophage polarization, immune subtype, molecular subtype, and immunomodulators, as well as the tumor immune microenvironment (TIME)-related prognosis, in LIHC. The study established six favorable prognostic models to predict LIHC prognosis based on the SOX4-associated genes. Finally, lenvatinib treatment can increase the expression of SOX4 in hepatocellular carcinoma cells and lead to drug resistance. Silencing SOX4 can effectively eliminate the drug resistance caused by lenvatinib treatment and inhibit the proliferation of cancer cells.Conclusions: This study highlights that SOX4 may serve as a promising therapeutic target for tumor treatment.


Introduction
Cancer ranks as a significant financial burden around the world.An estimated 19.3 million new cancer cases and almost 10.0 million cancer deaths occurred in 2020 worldwide [1].In the United States, the cases of global cancer burden are estimated to reach 28.4 million in 2040, a 47% increase from 2020 [1,2].In China, the age-standardized incidence rates (204.8 per 100,000) and age-standardized mortality rates (129.4 per 100,000) are higher than the global average level [3].In recent years, along with the development of gene expression profiling and large-scale cancer multi-omics databases, comprehensive analysis is available not only to generate a global view of different tumor types but also to provide insight into their biological characteristics and molecular mechanisms.This is crucial to facilitate tumor diagnosis and potential drug screening to aid in the development of precision medicine strategies of cancer.
Sex-determining region Y-related (SRY) high-mobility group (HMG) box 4 (SOX4) is a member of the C subgroup of SRY-related HMG box transcription factor family.It is comprised of three domains.The HMG-box domain (59-138 amino acids) can directly change the chromatin architecture and, subsequently, alter gene expression [4].The serine-rich region is a transactivation domain, whereas the glycine-rich region is a novel functional region [5].SOX4 diversely regulates many cellular biological events such as embryonic development and cell fate determination [6].SOX4 also plays an important role in tumorigenesis and cancer progression [7,8].In addition, SOX4 is closely related to the multidrug resistance of lung cancer, colon cancer, breast cancer, leukemia, and other malignant tumors [9][10][11][12].However, the specific drug resistance mechanism of SOX4 in LIHC is still unclear.These findings suggest SOX4 is a promising therapeutic target for cancer treatment.
In the present study, we firstly elucidated SOX4 expression and explored its correlation with prognosis in human tumor tissues.We also explored the genetic mutation and DNA promoter methylation level of SOX4 in different tumor types.Then, the correlation between the expression of SOX4 and immune cell infiltration, gene markers, and macrophage polarization across different tumor types, as well as immune score distribution, immune checkpoints-related gene, immunomodulator, immune, and molecular subtypes, were analyzed in liver hepatocellular carcinoma (LIHC).Moreover, we analyzed the association of SOX4 expression with DNA damage repair-, epithelial to mesenchymal transition (EMT)-, M6A methylation-, hypoxia-, ferroptosis-, and energy-metabolism-related genes in LIHC and constructed six prognostic nomograms to predict the 1-year, 3-year, and 5-year survival probabilities based on these SOX4-associated genes, respectively.Finally, the role of SOX4 in lenvatinib resistance and proliferation in hepatocellular carcinoma were identified in vitro.

Expression Analysis
TIMER2.0 (http://TIMER2.cistrome.org(accessed on 13 July 2021)) was used to study the expression of the SOX4 gene in tumor tissues and was matched to normal tissues.TIMER2.0 can analyze the tumor immunological, clinical, and genomic features across diverse cancer types."Gene_DE" module of TIMER2.0 can identify the up-regulated genes or downregulated genes in the tumors.Considering the limited availability of normal samples in the TCGA database, Gene Expression Profiling Interactive Analysis 2 (GEPIA 2) (http://GEPIA.cancer-pku.cn(accessed on 13 July 2021)) was used to study the differential gene expression of SOX4 between 9 tumors and adjacent normal tissues, including DLBC, LAML, LGG, TGCT, THYM, UCS, ACC, OV, and SARC.GEPIA 2 can integrate the TCGA database and GTEx database."Expression DIY-Box Plot" of GEPIA 2 can profile the tissuewise expression of one gene or a multigene signature using a box plot.The parameters were as follows: |Log2FC| cutoff = 1, p-value cutoff = 0.01, log scale = "Yes", jitter size = 0.04, and matched normal data = "Match TCGA normal and GTEx data"."Expression DIY-Stage Plot" of GEPIA 2 can profile the expression of one gene in different tumor stages.The parameters were as follows: use major stage = "Yes", and log scale = "Yes".The UALCAN cancer database (http://ualcan.path.uab.edu(accessed on 13 July 2021)) was used to study the differential total-protein SOX4 expression between 6 tumors and adjacent normal tissues, including breast cancer, ovarian cancer, colon cancer, clear cell renal cell carcinoma, uterine corpus endometrial carcinoma, lung adenocarcinoma, and pediatric brain cancer.Moreover, GEPIA 2 was also used to analyze the expression and survival role of lncRNA in LIHC.The "CPTAC analysis" module of UALCAN was used to study the protein expression using data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) Confirmatory/Discovery dataset.The Human protein atlas (HPA) (https://www.proteinatlas.org/(accessed on 13 July 2021)) was used to obtain immunohistochemistry images to map the protein expression levels of SOX4 in different tumor tissues.

Survival Analysis
GEPIA 2 and Kaplan-Meier Plotter (http://kmplot.com/analysis/(accessed on 13 July 2021)) were used to study the correlation of SOX4 expression with patient overall survival (OS) and disease-free survival (DFS) across different tumor types, respectively.Survival Analysis-Survival Plots of GEPIA 2 can perform survival analysis and plot a Kaplan-Meier curve.The parameters were as follows: group cutoff = median, cutoff-high (%) = 50, cutoff-low (%) = 50, hazard ratio (HR) = "Yes", 95% confidence interval = "Yes", and axis units = "Months".In this study, a "Kaplan-Meier plotter in pan-cancer" analysis was used to analyze the correlation between SOX4 expression and patient survival in 21 different cancers.The hazard ratio with 95% confidence intervals and log-rank p-value were also calculated and are shown in forest plots generated using GraphPad Prism 8.2.

Genetic Mutation Analysis
CBioPortal (https://www.cbioportal.org/(accessed on 13 July 2021)) were used to elucidate the potential role of SOX4 mutation in data extracted from TCGA dataset.The mutation frequency of SOX4 gene was analyzed with the "Cancer Types Summary" module of cBioPortal.The types, sites, case number and 3D structure of SOX4 mutation were analyzed with the "mutations" module of cBioPortal.The correlation of SOX4 mutation with patient OS, DFS, disease-specific survival (DSF), and progression-free survival (PFS) were analyzed with the "Comparison/Survival" module of cBioPortal.The correlation of SOX4 mutation with the fraction of the copy number of the altered genome was analyzed using the "Plots" module of cBioPortal.In addition, we analyzed the correlation of SOX4 expression with the tumor mutation burden (TMB) and microsatellite instability (MSI) across all tumors of TCGA using Spearman's correlation, and the correlation coefficient and p-value are shown in radar maps.

Methylation Analysis
The UALCAN cancer database was used to study the DNA promoter methylation level of SOX4 in tumor types.The "TCGA analysis" module of UALCAN was used to obtain the DNA promoter methylation level of SOX4.Especially, the correlation between the specific methylation site and SOX4 expression in LIHC was investigated with MEXPRESS (https://mexpress.be(accessed on 13 July 2021)).MEXPRESS can visualize the TCGA expression, DNA methylation, and clinical data, as well as the relationships among them.

Immune Characteristics Analysis
TIMER2.0 and TISIDB (http://cis.hku.hk/TISIDB/index.php(accessed on 13 July 2021)) were used to study the correlation of SOX4 expression with the immune characteristics of tumors.The "Gene" module of TIMER2.0 was used to study the correlation between SOX4 expression and the immune cells infiltration.The relationship between SOX4 expression and macrophage polarization was also studied.The "Gene_Corr" module of TIMER2.0 was used to study the correlation between SOX4 expression and the markers of immune infiltrating cells.In addition, the association between SOX4 expression and the subtypes immunomodulators was explored with the "Subtype" module of the TISIDB and the "Immunomodulator" module of the TISIDB in LIHC, respectively.

Prognostic Model Establishment
The raw counts of RNA-sequencing data (level 3) and corresponding clinical information of LIHC were obtained from the TCGA dataset.The least absolute shrinkage and selection operator (LASSO) regression technique was used for the predictor selection.Then, the multivariate analysis was used to identify the independent factors (two-sided p-value < 0.05) and generate a final nomograms.In addition, the LIHC patients were divided into two groups (high-risk and low-risk) based on the median risk score.The expression profiles of the prognostic genes in the two groups were visualized using a cluster heatmap.The predictive accuracy of the nomograms was measured with the area under receiver operating characteristic curve (AUC) of the receiver operating characteristic (ROC) curve.
Cancers 2023, 15, 5235 5 of 21 2.9.Cell Culture Human hepatocellular carcinoma cell lines HepG2 and Huh7 were purchased from Shanghai Cell Bank, Chinese Academy of Sciences.HepG2 cells were cultured in MEM medium containing 10% fetal bovine serum, supplemented with 2 mM glutamine, 100 U/mL penicillin and 100 µg/mL streptomycin.The culture medium of the Huh7 cells was DMEM containing 10% fetal bovine serum, and the final concentration was 100 U/mL penicillin and 100 µg/mL streptomycin.Both cells were cultured in a 5% CO 2 incubator at 37 • C.

Protein Extraction and Western Blotting
Total protein was extracted from the HepG2 and Huh7 cells using RIPA lysate containing protease inhibitors (Solarbio, Beijing, China).The cells were collected with a scraper and lysed on ice for 30 min.Centrifugation was conducted at 4 • C, 12,000 rpm for 30 min.After centrifugation, the supernatant was discarded and the protein concentration was quantified using the BCA protein concentration assay kit (Solarbio, Beijing, China).Protein samples were separated with SDS-PAGE and transferred to PVDF membrane (Millipore, Billerica, MA, USA).The membrane was blocked with 5% skimmed milk powder at 37 • C for 1 h.The cells were incubated overnight with anti-SOX4 (Affinity, DF2610, 1:1000) antibody.Finally, ImageJ software (National Institutes of Health, version 1.8.0.345) was used to analyze and process the protein strength.

Cell Proliferative Assay
After cell counting, an appropriate amount of fresh medium was added to prepare a cell suspension with a cell concentration of approximately 5 × 104/mL cells.The cells were seeded in 96-well plates, and 100 µL cell suspension was inoculated into each well, approximately 5000 cells per well.The 96-well plates were placed in an incubator at 37 • C and 5% CO 2 .After the cells were adherent, the cells were treated with lenvatinib or SOX4 silencing combined with lenvatinib (6 wells in each group) for 72 h.The old medium in the 96-well plate was discarded and replaced with a mixed solution of fresh medium and CCK8 (medium:CCK8 was 10:1).After incubation in an incubator for 2 h, the OD value of each well sample was measured at a wavelength of 450 nm with a multifunctional microplate reader, and the cell survival rate was calculated.

Trypan Blue Staining
Cells were seeded in culture flasks and placed in an incubator at 37 • C, 5% CO 2 .After the cells were adherent, the cells were treated with lenvatinib or knockout of SOX4 combined with lenvatinib and continued to incubate for 48 h.Trypsin was added to digest cells and prepare a single cell suspension (10 6 /mL cells).The cell suspension was mixed with 0.4% trypan blue solution at a ratio of 9:1.After three minutes, the living cells and dead cells were counted using a counting plate.

Detection of Apoptosis and Necrosis
Trypsin was added to the tumor cells treated with lenvatinib or SOX4 knockdown.After being centrifuged at 4 • C 1000 rpm for 5 min, the supernatant was discarded.The cells were resuspended in 1 mL of pre-cooled PBS, centrifuged at 1000 rpm at 4 • C for 5 min, and the supernatant was discarded.The cells were resuspended with 1 mL Cell Stain Buffer, followed by 5 µL Hoechst 33342 Stain and 5 µL PI Stain.Finally, the fluorescence intensity was detected using a fluorescence microscope after smear.

Statistical Analysis
R software version 4.0.0 and SPSS version 24.0 were used for the statistical analyses.The data are expressed as the means ± SD.The statistical analysis was performed using one-way ANOVA.A two-sided p-value < 0.05 was considered to be statistically significant.

Expression of SOX4 across Tumor Types
We first compared the SOX4 mRNA expression based on the TCGA data via the TIMER tool.As shown in Figure 1a, compared with the adjacent normal tissues, the expression of SOX4 was markedly increased in 18 tumor types or specific cancer subtypes including BLCA, BRCA, CHOL, COAD, ESCA, GBM, HNSC-HPV+, HNSC, KIRP, LIHC, LUAD, LUSC, PCPG, PRAD, READ, STAD, THCA, and UCEC, but it was significantly decreased in KICH and KIRC.There was no statistically significant difference in CESC, PAAD, and SKCM.Considering the limited availability of normal data in the TCGA database, we integrated the TCGA and GTEx data to study the difference in SOX4 mRNA expression in nine tumor types.As shown in Figure 1b, compared with the adjacent normal tissues, SOX4 expression was markedly increased in DLBC, LAML, LGG, TGCT, THYM, and UCS.However, there was no statistically significant difference in ACC, OV, and SARC.Moreover, we performed a protein expression analysis of SOX4 for breast cancer, colon cancer, ovarian cancer, clear cell renal cell carcinoma, UCEC, LUAD, and pediatric brain cancer using the "CPTAC analysis" module of UALCAN.As shown in Figure 1c, compared with the adjacent normal tissues, the protein expression of SOX4 was markedly increased in breast cancer, LUAD, and UCEC.There is no information available for SOX4 in other tumor types.The immunohistochemical staining of SOX4 in various tumors are shown in Supplementary Materials Figure S1.We also evaluated the correlation between the expression of SOX4 and the tumors pathological stages.SOX4 expression was found to significantly correlate with the pathological stages of ESCA, LIHC, PAAD, SKCM, and THCA (Figure 1d).

Prognostic Values of SOX4 across Tumor Types
We studied the prognostic values of SOX4 in patients.According to the GEPIA2 analysis results, high SOX4 expression was significantly correlated with poor OS in LGG and LIHC, and low SOX4 expression was significantly associated with poor OS in KIRC and THYM (Figure 2a).Moreover, high SOX4 expression was significantly correlated with poor DFS in ESCA and PAAD, and low SOX4 expression was significantly associated with poor DFS in UCEC (Figure 2b).According to the Kaplan-Meier analysis results, high SOX4 expression was significantly correlated with poor OS in LIHC, pancreatic ductal adenocarcinoma, and SARC.However, low SOX4 expression was significantly associated with poor OS in HNSC, ovarian cancer, and THYM (Supplementary Materials Figure S2a).For DFS, high SOX4 expression was significantly correlated with poor DFS in cervical squamous cell carcinoma, esophageal adenocarcinoma, LIHC, lung adenocarcinoma, pancreatic ductal adenocarcinoma, and THCA.However, low SOX4 expression was significantly associated with poor DFS in ovarian cancer and PCPG (Supplementary Materials Figure S2b).As shown in Figure 2c,d, the correlations of SOX4 expression with the OS and DFS of patients across 21 tumor types based on the Kaplan-Meier Plotter are exhibited in forest plots, respectively.

Genetic Mutation of SOX4 across Tumor Types
Genetic mutation is a major issue that contributes to tumorigenesis and prognosis.To elucidate the potential role of SOX4 mutation in tumors, we checked the mutation frequency of SOX4 gene across 32 human tumors using TCGA data from the cBioPortal database.As shown in Figure 3a, 24 of the 32 analyzed tumors had more than one mutation, and 8 tumors showed no mutation, including LAML, ACC, KICH, KIRP, MESO, PCPG, TGCT, and THCA.The patients with BLCA had the highest mutation frequency of the SOX4 gene (16.79%),followed by OV (6.16%).Among all types of genetic mutations, the "amplification" type was the most frequent alteration.Notably, the "amplification" type was the sole alteration of the SOX4 gene in the patients with CHOL, DLBC, LIHC, OV, PAAD, and UVM.In addition, the "mutated" type was the sole alteration of the SOX4 gene in the patients with KIRC and THYM.The types, sites, and case number of the SOX4 mutation (including missense, truncating, and in-frame) are presented in Figure 3b.The frequency of the somatic mutation was 0.3%.There were 35 mutations sites between 0 and 474 amino acids.Missense was the main type of genetic mutation.As shown in Figure 3c, compared with the samples without any alteration of the SOX4 gene (i.e., unaltered group, number of cases = 5277), the samples with at least one alteration of the SOX4 (i.e., altered group, number of cases = 106) showed poorer prognosis with regard to DFS but not PFS, OS, and DSS.As shown in Supplementary Materials Figure S3, the mutation count of the SOX4 gene was positively correlated with the fraction of the copy number altered genome (Spearman = 0.34, Pearson = 0.21).Meanwhile, we studied the correlation of SOX4 expression with TMB and MSI.The results indicate that the aberrant expression of SOX4 was positively associated with TMB in the patients with BLCA, LUAD, PAAD, PRAD, and THCA, but it was negatively correlated with that of COAD and THYM (Figure 3d, Supplementary Materials Figure S4).The aberrant expression of SOX4 was positively associated with MSI in the patients with LUSC, READ, and UCEC and negatively correlated with that of COAD, PRAD, and SKCM (Figure 3e, Supplementary Materials Figure S5).

Prognostic Values of SOX4 across Tumor Types
We studied the prognostic values of SOX4 in patients.According to the GEPIA2 analysis results, high SOX4 expression was significantly correlated with poor OS in LGG and LIHC, and low SOX4 expression was significantly associated with poor OS in KIRC and THYM (Figure 2a).Moreover, high SOX4 expression was significantly correlated with poor DFS in ESCA and PAAD, and low SOX4 expression was significantly associated with poor DFS in UCEC (Figure 2b).According to the Kaplan-Meier analysis results, high squamous cell carcinoma, esophageal adenocarcinoma, LIHC, lung adenocarcinoma, pan creatic ductal adenocarcinoma, and THCA.However, low SOX4 expression was signif cantly associated with poor DFS in ovarian cancer and PCPG (Supplementary Materia Figure S2b).As shown in Figure 2c,d

ceRNA Regulatory Network of SOX4 in LIHC
Increasing evidence shows that changes to and dysfunction of lncRNA lead to abnormal gene expression and promote the formation, progression, and metastasis of many types of cancer.lncRNAs can absorb miRNA and then promote mRNA expression.In order to explore the relationship between SOX4 and lncrna and miRNA we, therefore, analyzed the ceRNA regulatory network of SOX4 in tumor tissues of LIHC.We firstly predicted the potential miRNA of SOX4 with six target prediction databases, including DIANA, PITA, TargetScan, miRTarBase, miRmap, and mirDIP.As shown in the Venn diagram (Figure 4a), there are 215 overlapping miRNAs identified in at least four of the six target prediction databases.Among these miRNAs, only two miRNAs (hsa-miR-139-3p and hsa-miR-30e-5p) were significantly lower than that of the adjacent normal tissues and negatively correlated with the expression of SOX4 in LIHC (Figure 4b,c).Moreover, has-miR-139'3p's expression was significantly correlated with poor OS in LIHC.However, there was no statistically significant difference between the expression of hsa-miR-30e-5p and OS in LIHC (Figure 4d).Therefore, we further focused on the role of hsa-miR-139-3p in LIHC and predicted its upstream lncRNAs.There were 476 lncRNAs identified using LncBase Predicted v.2 database.Among these lncRNAs, only LINC00152's expression was significantly higher than that of the adjacent normal tissue and was positively correlated with the expression of SOX4 in LIHC (Figure 4e,f).Moreover, LINC00152's expression was correlated with poor OS in LIHC (Figure 4g).LINC00152's expression was also negatively correlated with the expression of hsa-miR-139-3p in LIHC.So far, we found that LINC00152 could function as a ceRNA to regulate SOX4 expression by sponging hsa-miR-139-3p in LIHC.

Immune Characteristics of SOX4 in LIHC
The tumor immune microenvironment (TIME) is involved in tumor clonal evolution, growth, metastasis, prognosis, and drug resistance, as well as therapeutic outcome.Therefore, we further analyzed the immune characteristics of SOX4 in LIHC.First, we assessed the relationship between SOX4 expression and the infiltration level of immune cells.The

Immune Characteristics of SOX4 in LIHC
The tumor immune microenvironment (TIME) is involved in tumor clonal evolution, growth, metastasis, prognosis, and drug resistance, as well as therapeutic outcome.Therefore, we further analyzed the immune characteristics of SOX4 in LIHC.First, we assessed the relationship between SOX4 expression and the infiltration level of immune cells.The correlations between SOX4 expression and the immune infiltrating cells (including CD8+ T cells, CD4+ T cells, B cells, neutrophils, macrophages, and myeloid dendritic cells) across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S6).Especially in LIHC, the result shows that SOX4 expression was significantly correlated with tumor purity in TIMER.Notably, SOX4 expression was positively correlated with CD4+ T cells, B cells, neutrophils, macrophages, and myeloid dendritic cells but not CD8+ T cells (Figure 5a).For the macrophage polarization analysis, SOX4 expression was positively correlated with M0 polarization but not with M1 or M2 polarization (Figure 5b).Second, we analyzed the relationship between SOX4 expression and the markers of CD8+ T cells, B cells, neutrophils, macrophages, myeloid dendritic cells, tumor-associated macrophages, monocytes, and natural killer cells, as well as Tfh, Th1, Th2, Th9, Th17, Th22, Treg, and exhausted T cells.The correlations between SOX4 expression and the markers of the above cells across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S7).Especially in LIHC, the expression of SOX4 was positively associated with CD19, CD38, CD8A, CD8B, CXCR5, ICOS, BCL6, CCR3, GATA3, TGFBR2, IRF4, IL21R, IL23R, STAT3, CCR10, AHR, CCR8, CTLA4, CD68, CD80, CD86, XCL1, CD7, MPO, and CD1C and negatively associated with ARG1 and CD14.

Prognostic Models Based on SOX4-Associated Genes in LIHC
We further established six multigene prognostic models to predict LIHC prognosis based on SOX4-associated genes, respectively.The genes were initially collected through a comprehensive literature search.For the prognostic model based on the SOX4-associated DNA damage repair-related genes, 200 genes were collected and their correlations with SOX4 expression across various tumor types were visualized using cluster heatmaps

Prognostic Models Based on SOX4-Associated Genes in LIHC
We further established six multigene prognostic models to predict LIHC prognosis based on SOX4-associated genes, respectively.The genes were initially collected through a comprehensive literature search.For the prognostic model based on the SOX4-associated DNA damage repair-related genes, 200 genes were collected and their correlations with SOX4 expression across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S14a).In LIHC, 21 genes had nonzero coefficients in the LASSO analysis (Supplementary Materials Figure S14b,c).Multivariate logistic regression indicated that RFC3, RAD54B, MUTYH, MGMT, HAP, and UVSSA were independent risk factors for OS (Table 1).These genes were used to construct the nomogram (Figure 6a), where risk score = (−0.1752)× RFC3 + (1.4157) × RAD54B + (0.3664) × MUTYH + (−0.3797) × MGMT + (0.406) × HAP1 + (−0.3656) × UVSSA.The risk scores of the high-risk group and low-risk group are shown in Supplementary Materials Figure S14d.The survival time and survival status of the patients are shown in Supplementary Materials Figure S14e.The expression profiles of the prognostic genes are visualized in Supplementary Materials

Role of SOX4 Knockdown in Lenvatinib-Treated LIHC Cells
To further validate the results of SOX4 in above pan-cancer analysis and LIHC, analyzed the role of SOX4 knockdown in lenvatinib-treated LIHC cells in vitro.As kno For the prognostic model based on the SOX4-associated EMT-related genes, 95 genes were collected and their correlations with SOX4 expression across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S15a).In LIHC, nine genes had nonzero coefficients in the LASSO analysis (Supplementary Materials Figure S15b,c).A multivariate logistic regression indicated that ACTA2, PHLDA2, SGCB, and NKX3-2 were independent risk factors for OS (Table 1).These genes were used to construct the nomogram (Figure 6b), in which the risk score = (−0.2898)× ACTA2 + (0.2364) × PHLDA2 + (0.3189) × SGCB + (0.5039) × NKX3-2.The risk scores, survival time and survival status, expression profiles, Kaplan-Meier survival analysis, and ROC curve analyses are shown in Supplementary Materials Figure S15d-h.A multivariate logistic regression indicated that ZC3H13 and YTHDF2 were independent risk factors for OS (Table 1).These genes were used to construct the nomogram (Figure 6c), in which risk score = (−0.3353)× ZC3H13 + (1.0005) × YTHDF2.
For the prognostic model based on the SOX4-associated hypoxia-related genes, 75 genes were collected and their correlations with SOX4 expression across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S16a).In LIHC, 12 genes had nonzero coefficients in the LASSO analysis (Supplementary Materials Figure S16b,c).A multivariate logistic regression indicated that CUL2, EPO, and UBB were independent risk factors for OS (Table 1).These genes were used to construct the nomogram (Figure 6d), in which risk score = (0.6505) × CUL2 + (0.1865) × EPO + (−0.2321) × UBB.The risk scores, survival time and survival status, expression profiles, Kaplan-Meier survival analysis, and ROC curve analyses are shown in Supplementary Materials Figure S16d-h.
For the prognostic model based on the SOX4-associated energy-metabolism-related genes, 100 genes were collected and their correlations with SOX4 expression across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S17a).In LIHC, 10 genes had nonzero coefficients in the LASSO analysis (Supplementary Materials Figure S17b,c).A multivariate logistic regression indicated that GGT3P, LDHA, and NQO2 were independent risk factors for OS (Table 1).These genes were used to construct the nomogram (Figure 6e), in which risk score = (2.395)× GGT3P + (0.5369) × LDHA + (0.2521) × NQO2.The risk scores, survival time and survival status, expression profiles, Kaplan-Meier survival analysis, and ROC curve analyses are shown in Supplementary Materials Figure S17d-h.
For the prognostic model based on the SOX4-associated ferroptosis-related genes, 24 genes were collected and their correlations with SOX4 expression across various tumor types were visualized using cluster heatmaps (Supplementary Materials Figure S18a).In LIHC, nine genes had nonzero coefficients in the LASSO analysis (Supplementary Materials Figure S18b,c).Multivariate logistic regression indicated that SAT1, SLC7A11, and CISD1 were independent risk factors for OS (Table 1).These genes were used to construct the nomogram (Figure 6f), in which risk score = (−0.2718)× SAT1 + (0.2648) × SLC7A11 + (0.4011) × CISD1.The risk scores, survival time and survival status, expression profiles, Kaplan-Meier survival analysis, and ROC curve analyses are shown in Supplementary Materials Figure S18d-h.

Role of SOX4 Knockdown in Lenvatinib-Treated LIHC Cells
To further validate the results of SOX4 in above pan-cancer analysis and LIHC, we analyzed the role of SOX4 knockdown in lenvatinib-treated LIHC cells in vitro.As known, lenvatinib was approved as a first-line drug for the treatment of unresectable hepatocellular carcinoma [13].Firstly, Huh7 and HepG2 cells were treated with 0, 5, or 10 uM lenvatinib for 24 or 48 h, respectively.The results show that lenvatinib increased the expression of SOX4 in HepG2 and Huh7 cells in a time-and dose-dependent manner (Figure 7a).In order to further explore the relationship between SOX4 expression and lenvatinib resistance, we knocked down SOX4 in LIHC cells.HepG2 or Huh7 cells treated with lenvatinib (10 uM) were treated with SOX4 siRNA for 48 h to analyze cell proliferation and necrosis.The CCK8 assay and trypan blue staining showed that, compared with the lenvatinib treatment group, lenvatinib combined with SOX4 silencing significantly inhibited the proliferation and viability of HepG2 and Huh7 cells (Figure 7b,c).The apoptosis and necrosis experiments also showed that, compared with lenvatinib treatment alone, lenvatinib combined with SOX4 silencing significantly increased the necrosis of HepG2 and Huh7 cells (Figure 7d).
Cancers 2023, 15, x FOR PEER REVIEW 18 of 23 lenvatinib treatment group, lenvatinib combined with SOX4 silencing significantly inhibited the proliferation and viability of HepG2 and Huh7 cells (Figure 7b,c).The apoptosis and necrosis experiments also showed that, compared with lenvatinib treatment alone, lenvatinib combined with SOX4 silencing significantly increased the necrosis of HepG2 and Huh7 cells (Figure 7d).

Discussion
A comprehensive analysis is of great importance for the accurate diagnosis and effective therapy of tumors [14,15].This is crucial not only to tackle the biological characteristics and molecular mechanisms of different cancer types but also to provide some insight for researchers into the development of more precise therapeutic strategies against various cancers.Therefore, many researchers are focusing on exploring ways to generate a global view of different cancer types [16,17].The present study aimed to visualize the expression, survival, mutations, methylation, ceRNA network, immunity, and prognostic models, as well as explore the potential role, of SOX4 across different tumor types and specific LIHC using different online tools.
SOX4 is a critical transcription factor that is involved in many cellular events, such as stemness, differentiation, progenitor development, and tumorigenesis, contributing to diverse pathological conditions [4,18].Although emerging research has evaluated the expression level and function of SOX4 in certain tumors, its roles across various tumor types, especially aspects related to prognostic potential and clinical significance, have not been systematically studied.In this study, we comprehensively analyzed the expression levels of SOX4 between tumors and normal tissues for the first time.Our study indicated that the mRNA expression of SOX4 was upregulated in most of the tumors.In addition, the study indicated that the protein expression of SOX4 was significantly higher in breast cancer, LUAD, and UCEC.It was also found that the expression of SOX4 was correlated with the pathological stages of ESCA, LIHC, PAAD, SKCM, and THCA.Overall, these results strongly suggest that SOX4 might exert specific and even contrasting functions in different tumor types.
To further clarify the functions of SOX4 in different tumor types, this study analyzed the correlation between SOX4 expression and patients' prognosis.Our study shows that SOX4 upregulation was associated with a poor prognostic outcome in several tumors, especially in LIHC.Therefore, SOX4 might be a potential prognostic and diagnostic marker of survival outcomes.Indeed, the prognostic values of SOX4 have previously been reported in some malignant tumors.Some experimental studies have revealed that SXO4 expression is correlated with poor prognosis in patients with tumors, such as CHOL [19] and LAML [20].The prognostic value of SOX4 was partly validated again in our study, despite some controversial results.The role of SOX4 in these cancers is still open to question and further discussion.
To clarify the underlying mechanisms that drive the observed results, we further checked the genetic mutations and DNA promoter methylation level of SOX4.At the gene level, we found that the samples with at least one alteration of the SOX4 gene showed poor prognosis.In addition, the expression of SOX4 was significantly associated with MSI and TMB status.At the DNA level, the promoter methylation level of SOX4 was upregulated in ESCA, KIRC, KIRP, LIHC, LUSC, PAAD, SARC, and TGCT, while it was downregulated in BLCA, BRCA, HNSC, THCA, and UCEC.Especially in LIHC, a strong correlation between SOX4 expression and the methylation sites was observed.Indeed, the accumulation of genetic mutations and epigenetic alterations play a pivotal role during tumorigenesis.It is well known that many tumors harbor several genetic mutations.Genetic mutations drive the development of tumors and affect tumor prognosis [21].Epigenetic alterations regulate gene expression and establish cell-type-specific temporal and spatial expression patterns [22].These may partially explain the complex roles of SOX4 in the prognosis of various tumors.ceRNA regulatory networks are recognized as important regulators of gene expression at the posttranscriptional level.Accumulating evidence indicates that ceRNA regulatory networks regulate many biological processes, especially in tumors [23].In our study, we firstly predicted the potential miRNAs of SOX4.We found that hsa-miR-139-3p was not only significantly decreased but also negatively correlated with SOX4 expression and poor OS in LIHC.The previously reported miR-363-3P and miR-129-2 were also found to be downregulated in HCC and negatively regulate SOX4 levels in vitro [24,25].Indeed, hsa-miR-139-3p functions as a biomarker in different cancer types, such as colorectal cancer and hepatocellular carcinoma [26].We further predicted the upstream lncRNAs of hsa-miR-139-3p.We found that LINC00152 was not only significantly higher than that of the adjacent normal tissues but also positively correlated with SOX4 expression and poor OS in LIHC.Accumulated evidence shows that LINC00152 plays an important role in carcinogenesis by disturbing various signaling pathways in several cancers [27].Nevertheless, no previous study has evaluated the importance of the role of the LINC00152/hsa-miR-139-3p/SOX4 ceRNA regulatory network in LIHC.In our study, we evaluated the prognostic values of LINC00152/hsa-miR-139-3p/SOX4, which will provide novel insight into the treatment of LIHC.
The TIME was another crucial aspect of this study.The TIME acted as a crucial regulator in the development, progression, and immune escape of various types of cancer [28].A better definition and understanding of the TME could open novel avenues for the curative treatment of metastatic tumors.In this study, the analysis revealed that SOX4 expression was associated with CD4+ T cells, B cells, neutrophils, macrophages, and myeloid dendritic cells but not CD8+ T cells.Moreover, SOX4 expression was correlated with the M0 polarization of macrophages.We also confirmed the relationship between SOX4 expression and the markers of immune-infiltrating cells, especially in LIHC.We further investigated the associations between SOX4 expression and the immune subtypes, molecular subtypes, and immunomodulators in LIHC.We further established a multigene prognostic model to predict the clinical prognosis of LIHC.Notably, the model also displayed excellent predictive discrimination.Thus, this model is a sensitive prediction tool under the promise of guaranteeing accuracy.
It is well known that tumors are shaped by combined action of lifestyle, environmental, eating habits, physical activity, and genetic factors.Some reliable models that can predict the prognosis are urgently required.As known, aberrant DNA damage, EMT, M6A methylation, hypoxia, energy metabolism, and ferroptosis are important events generally involved in tumor onset and development [29][30][31][32][33][34].Therefore, we established six multigene prognostic models to predict LIHC prognosis based on the SOX4-associated genes, including DNA damage repair-related genes, EMT-related genes, M6A methylation-related genes, hypoxia-related genes, energy-metabolism-related genes, and ferroptosis-related genes.The estimated 1-, 3-, and 5-year survival probabilities could easily be calculated using the nomograms.All models also displayed excellent predictive discrimination with regard to the AUCs for the 1-, 3-, and 5-year survival.To date, there exists no study that has evaluated the prognostic values of these genes in patients with LIHC.The present study might be the first to have constructed nomograms for the prediction of LIHC prognosis based on the systematic assessment of cellular core genes.
Lenvatinib is a multitarget tyrosine kinase inhibitor that inhibits the growth and angiogenesis of malignant cells by inhibiting the activation of multiple receptor tyrosine kinase signaling pathways, thereby effectively treating a variety of cancers.In 2018, lenvatinib was approved for first-line treatment in patients with advanced hepatocellular carcinoma in the United States, the European Union, Japan, and China [35].However, some patients developed drug resistance after using lenvatinib for a period of time.Increasing evidence shows that SOX4 is closely related to tumor drug resistance.In our study, we found that lenvatinib treatment upregulated the expression of SOX4 in Huh7 and HepG2 cells in a dose-and time-dependent manner, indicating that SOX4 may associated with lenvatinib drug resistance in LIHC.In addition, we also demonstrated that, compared with lenvatinib treatment, lenvatinib treatment combined with SOX4 silencing significantly inhibited the proliferation of LIHC cells and increased necrosis.

Conclusions
In conclusion, SOX4 was widely overexpressed in tumor tissues and associated with unfavorable prognoses, genetic mutation, and DNA methylation level, especially in LIHC.Moreover, the results provide novel evidence supporting LINC00152/hsa-miR-139-3p/SOX4 as a crucial target for the treatment of LIHC.The results also provide insight into the significant role of SOX4 expression in immune cell infiltration, macrophage polarization, immune subtype, molecular subtype, and immunomodulators, as well as TIME-related prognosis, in LIHC.Furthermore, this study established six favorable prognostic models to predict LIHC prognosis based on the SOX4-associated genes, including DNA damage repair-related genes, EMT-related genes, M6A methylation-related genes, hypoxia-related genes, energy-metabolism-related genes, and ferroptosis-related genes.Finally, we found that SOX4 played an important role in the drug resistance of lenvatinib in LICH in vitro.Altogether, this study emphasizes the critical roles of SOX4 in the diagnosis and prognosis of tumors, especially in LIHC, and as a promising therapeutic target for tumor treatment.Informed Consent Statement: All data for this study are available from the TCGA and GEO databases.The TCGA and GEO belong to public databases.The patients involved in the database obtained ethical approval.Users can download relevant data for free for research and publishing relevant articles.Our study was based on open source data, so there are no ethical issues and other conflicts of interest.

Figure 1 .
Figure 1.Expression of SOX4 across tumor types: (a) expression of SOX4 gene in different tumors or specific tumor subtypes determined using TIMER2; (b) expression of SOX4 gene in different tumors determined using GEPIA 2; (c) protein expression of SOX4 gene in different tumors determined using UALCAN; (d) correlation between SOX4 expression and the different pathological stages of tumors determined using GEPIA2.* p < 0.05, ** p < 0.01, and *** p < 0.001.

Figure 1 .
Figure 1.Expression of SOX4 across tumor types: (a) expression of SOX4 gene in different tumors or specific tumor subtypes determined using TIMER2; (b) expression of SOX4 gene in different tumors determined using GEPIA 2; (c) protein expression of SOX4 gene in different tumors determined using UALCAN; (d) correlation between SOX4 expression and the different pathological stages of tumors determined using GEPIA2.* p < 0.05, ** p < 0.01, and *** p < 0.001.
, the correlations of SOX4 expression with the OS an DFS of patients across 21 tumor types based on the Kaplan-Meier Plotter are exhibited i forest plots, respectively.

Figure 2 .
Figure 2. Prognostic values of SOX4 across tumor types: (a) correlation of SOX4 expression with O of different tumors in the TCGA determined using GEPIA2; (b) correlation of SOX4 expression wit DFS of different tumors in the TCGA determined using GEPIA2; (c) forest plot of the correlation o SOX4 expression with OS across 21 tumor types determined using Kaplan-Meier Plotter; (d) fore plot of the correlation of SOX4 expression with DFS across 21 tumor types determined by Kaplan Meier Plotter.

Figure 2 . 23 Figure 3 .Figure 3 .
Figure 2. Prognostic values of SOX4 across tumor types: (a) correlation of SOX4 expression with OS of different tumors in the TCGA determined using GEPIA2; (b) correlation of SOX4 expression with DFS of different tumors in the TCGA determined using GEPIA2; (c) forest plot of the correlation of SOX4 expression with OS across 21 tumor types determined using Kaplan-Meier Plotter; (d) forest plot of the correlation of SOX4 expression with DFS across 21 tumor types determined by Kaplan-Meier Plotter.

Cancers 2023 , 23 Figure 5 .
Figure 5.Immune characteristics of SOX4 in LIHC: (a) correlation of SOX4 expression with immune cell infiltration; (b) correlation of SOX4 expression with macrophage polarization; (c) correlation of SOX4 expression with immune subtype; (d) correlation of SOX4 expression with molecular subtype; (e) nomogram constructed based on ten SOX4-associated immune-related genes to predict LIHC prognosis.

Figure 5 .
Figure 5.Immune characteristics of SOX4 in LIHC: (a) correlation of SOX4 expression with immune cell infiltration; (b) correlation of SOX4 expression with macrophage polarization; (c) correlation of SOX4 expression with immune subtype; (d) correlation of SOX4 expression with molecular subtype; (e) nomogram constructed based on ten SOX4-associated immune-related genes to predict LIHC prognosis.

Figure S14f .Figure 6 .
Figure S14f.The Kaplan-Meier survival analysis is shown in Supplementary Materials Figure S14g.The ROC curve analyses are shown in Supplementary Materials Figure S14h.cers 2023, 15, x FOR PEER REVIEW 17 o

Figure 6 .
Figure 6.Prognostic models based on SOX4-associated genes in LIHC: (a) nomogram constructed based on six SOX4-associated DNA damage repair-related genes to predict LIHC prognosis; (b) nomogram constructed based on four SOX4-associated EMT-related genes to predict LIHC prognosis; (c) nomogram constructed based on two SOX4-associated M6A methylation-related genes to predict LIHC prognosis; (d) nomogram constructed based on three SOX4-associated hypoxia-related genes to predict LIHC prognosis; (e) nomogram constructed based on three SOX4-associated energymetabolism-related genes to predict LIHC prognosis; (f) nomogram constructed based on three SOX4-associated ferroptosis-related genes to predict LIHC prognosis.

Figure 7 .
Figure 7. SOX4 silencing increased the sensitivity of lenvatinib to hepatocellular carcinoma cells and cell necrosis: (a) Western blotting of SOX4 in Huh7 or HepG2 cells treated with 0, 5, and 10 uM lenvatinib for 24 or 48 h; (b) cell viability of Huh7 and HepG2 cells treated with 10 uM lenvatinib and lenvatinib combined with knockdown of SOX4 was detected; (c) trypan blue staining of Huh7 and HepG2 cells treated with lenvatinib combined with knockdown of SOX4; (d) apoptosis and necrosis staining of Huh7 and HepG2 cells after treatment with lenvatinib combined with knockdown of SOX4 for 48 h.* p < 0.05, ** p < 0.01, *** p < 0.001.

Figure 7 .
Figure 7. SOX4 silencing increased the sensitivity of lenvatinib to hepatocellular carcinoma cells and cell necrosis: (a) Western blotting of SOX4 in Huh7 or HepG2 cells treated with 0, 5, and 10 uM lenvatinib for 24 or 48 h; (b) cell viability of Huh7 and HepG2 cells treated with 10 uM lenvatinib and lenvatinib combined with knockdown of SOX4 was detected; (c) trypan blue staining of Huh7 and HepG2 cells treated with lenvatinib combined with knockdown of SOX4; (d) apoptosis and necrosis staining of Huh7 and HepG2 cells after treatment with lenvatinib combined with knockdown of SOX4 for 48 h.* p < 0.05, ** p < 0.01, *** p < 0.001.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10.3390/cancers15215235/s1, FigureS1: Representative immunohistochemical staining results of SOX4 in various types of tumors, Figure S2: Prognostic values of SOX4 determined by Kaplan-Meier Plotter, Figure S3: The correlation of SOX4 mutation count with the fraction of copy number altered genome, Figure S4: The correlation of SOX4 expression with TMB in different tumors, Figure S5: The correlation of SOX4 expression with MSI in different tumors, Figure S6: The correlation between SOX4 expression and the infiltration level of immune cells visualized by the cluster heatmaps, Figure S7: The correlation between SOX4 expression and the markers of immune and immune-related cells visualized by the cluster heatmap, Figure S8: The correlation between SOX4 expression and immune subtypes, molecular subtypes in different tumors, Figure S9: The correlation between SOX4 expression and three kinds of immunomodulators in different tumors visualized by the cluster heatmap, Figure S10: The correlation of SOX4 expression with immunoinhibitors visualized by the box plots, Figure S11: The correlation of SOX4 expression with immunostimulators visualized by the box plots, Figure S12: The correlation of SOX4 expression with MHC molecules visualized by the box plots, Figure S13: The multigene prognostic model based on SOX4 associated Immunerelated genes to predict LIHC prognosis, Figure S14: The multigene prognostic model based on SOX4 associated DNA damage-related genes to predict LIHC prognosis, Figure S15: The multigene prognostic model based on SOX4 associated EMT-related genes to predict LIHC prognosis, Figure S16: The multigene prognostic model based on SOX4 associated hypoxia-related genes to predict LIHC prognosis, Figure S17: The multigene prognostic model based on SOX4 associated energy metabolismrelated genes to predict LIHC prognosis, Figure S18: The multigene prognostic model based on SOX4 associated ferroptosis-related genes to predict LIHC prognosis.Author Contributions: Methodology, X.D.; software, Y.W.; validation, H.G.; resources, Q.W.; writingoriginal draft preparation, X.D.; writing-review and editing, H.W.; visualization, Y.W.; supervision, S.R.All authors have read and agreed to the published version of the manuscript.Funding: This research was supported by the Natural Science Foundation of Hebei province (No.H2020307019), 2022 government-funded clinical talents project of Hebei Province (No. 361003), and 333 Talent Project of Hebei Province (No. A202102012).