HS3ST3A1 and CAPN8 Serve as Immune-Related Biomarkers for Predicting the Prognosis in Thyroid Cancer

Background Thyroid cancer (TC) tends to be a common malignancy worldwide and results in various outcomes due to its different subtypes. The tumor microenvironment (TME) was demonstrated to play crucial roles in various malignancies, including thyroid cancer. This study combined the ESTIMATE and CIBERSORT algorithms, identified four TME-related genes, and evaluated their correlation with clinical characteristics. These findings revealed the malignant performance of TME in TC, and the TME-related DEGs might serve as prognostic biomarkers, which can be utilized for the prediction of immunotherapy effects in patients with TC. Methods The clinical and gene expression profiles of TC patients were collected from the TCGA dataset. The ESTIMATE algorithm was utilized to estimate stromal and immune scores and predict the level of stromal and immune cell infiltration. The differential expressed genes related to TME were filtered by the “limma” package in R, and the PPI network was constructed by a string website. KEGG pathway and GO analyses were performed to investigate the biological progression and molecular functions of TME-related DEGs. Then, univariate Cox regression analysis was employed to screen four genes correlated with clinical characteristics. GSEA was conducted to assess their roles in the TME of TC. To further investigate the association between TME-related genes and tumor-infiltrating immune cells (TIICs), the CIBERSORT algorithm was performed. Finally, the malignancy behaviors of the two genes were verified by RT-qPCR, IHC, MTT, colony formation, and transwell assays. Results Four TME-related DEGs, LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8, were identified and were significantly predictive of poor overall survival. KEGG and GO pathway analysis established that the TME-related DEGs were involved in immune responses and pathways in cancer. Furthermore, the malignancy behaviors of HS3ST3A1 and CAPN8 were verified by cellular functional experiments. These results revealed that the TME-related genes HS3ST3A1 and CAPN8 were able to serve as predictors of prognosis in patients with TC. Conclusion HS3ST3A1 and CAPN8 may serve as valuable prognostic biomarkers and TME indicators, which can be utilized for the prediction of immunotherapy effects and provide novel treatment strategies for patients with TC.


Introduction
Tyroid cancer (TC) occurs commonly worldwide and represents 3% of the global incidence of all cancers, which shows a continuously increasing in the past three decades with 586000 new patients in 2020 [1,2]. In females, it develops three or four times more frequently but is less destructive than in males [3,4]. Te classifcation of TC subtypes is according to their diferentiated degree, which was commonly classifed into 4 histological types. Follicular thyroid carcinoma (FTC) and papillary thyroid carcinoma (PTC) are well-diferentiated thyroid cancers (WDTC) with better prognosis, while there are still poorly-diferentiated thyroid cancer (PDTC) and anaplastic thyroid carcinoma (ATC) [5]. PTC is the most common subtype that occupies nearly 80% of DTC [6]. Te initial management of DTC includes surgical resection, radioactive iodine ablation (RAI), and thyroid-stimulating hormone (TSH) suppression, after which, patients usually owe good outcomes. However, 10-15% of patients with thyroid cancer have recurrent disease, 5% develop distant metastasis (lungs and bones), and cancer-specifc death occurs in some cases [7]. Tus, it is essential to explore emerging biomarkers that contribute to prognosis prediction in thyroid cancer.
In previous studies, it has been demonstrated that the tumor microenvironment (TME) is intimately correlated to the development and prognosis of multicancer types [8][9][10]. Terefore, increased attention has been paid to the role of initiation and progression that TME plays in cancer. Stomal and immune cells stand for two primary nontumor members in TME, while there are still other complex components, including extracellular matrix (ECM) and infammatory mediators [11]. Emerging evidence has indicated that tumor-infltrating immune cells (TIICs), for instance, regulatory T cells (Tregs), CD8+ T cells, CD4+ T cells, and tumor-associated macrophages (TAMs), are involved in various malignancies [12,13]. Despite the previous study identifying several critical genes related to TME in thyroid cancer, further investigations are still necessary [14].
Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) is a standard algorithm for calculating scores of immune and stromal. Here, we utilized the ESTIMATE algorithm to calculate immune and stromal scores for patients with TC based on the gene expression profles from Te Cancer Genome Atlas (TCGA). Ten, the core genes involved in clinical outcomes were identifed via univariate Cox regression analysis. In this way, we expect to fnd out emerging genes, which can perform as biomarkers of prognostic prediction for patients with TC.

Data Collection and Mining.
Te gene expression and clinical profles of TC patients were collected from the TCGA dataset (https://portal.gdc.cancer.gov/). Te fltration conditions were set as "thyroid gland" and "TCGA-THCA". Finally, 568 patients with TC (58 normal specimens and 510 cancer specimens) were selected. Te composition of the TME in these specimens was evaluated by the ESTIMATE algorithm, and the results were shown as stromal score, immune score, and estimate score. "limma" package in R was performed to screen diferential expressed genes (DEGs) with criteria as follow: (1) |log 2FC| > 1; (2) FDR <0.05. After that, the "heatmap" package in R was employed to make TME-related DEGs visualization as heatmaps, and the "VennDiagram" package in R was employed to select similar genes in immune and stromal cells.

Functional Analysis of DEGs. Kyoto Encyclopedia of
Genes and Genomes (KEGG) and Gene Ontology (GO) analysis had been employed in TME-related genes. Te "clusterProfler," "enrichplot," and "ggplot2" packages were utilized for the visualization of KEGG and GO analysis results, which showed the biological process, molecular function, and pathways of TME-related DEGs. Data with P < 0.05 and q < 0.05 were considered with statistical signifcance.

PPI Network and Cox
Analysis of DEGs. Te STRING database (https://cn.string-db.org/) was employed to construct an interaction network of TME-related DEGs. Interactions with integrated scores higher than 0.95 were selected and visualized by Cytoscape (version 3.8.2). Te "survival" package in R was utilized to perform univariate Cox regression to determine the association of DEGs with the prognosis of TC. Ultimately, four core genes were screened to perform the following analysis.

Analysis of Tumor Immunoreaction.
To investigate the capable pathway that TME-related DEGs act on in the TME of TC, the Gene Set Enrichment Analysis (GSEA) was conducted using GSEA version 4.1.0 (Broad Institute, Cambridge, MA, United States). Te statistical signifcance was set as NOM P value <0.05 and FDR <0.25. Te correlation between TME-related DEGs and TIICs was assessed by the CIBERSORT algorithm, which is a deconvolution algorithm based on RNA-seq data to estimate the proportion of 22 immune cells in each specimen.

Correlation between TME-Related DEGs and
Clinicopathological Characteristics. Te overall survival (OS) was compared between the TC samples with high/low expression of TME-related DEGs through Kaplan-Meier survival analysis. Furthermore, the correlation analysis between DEGs and clinicopathological characteristics was performed, and results with P < 0.05 and q < 0.05 were considered with statistical signifcance.

Transfection of Short Hairpin RNAs (shRNAs).
Te shRNAs of HS3ST3A1 and CAPN8, as well as their control groups, were obtained from RiboBio (Guangzhou, China), which sequences were exhibited in the supplementary fle:

MTTand Colony Formation
Assay. Te MTT assay was performed by seeding 2 * 10 3 cells in 96-well plates after transfection for 48 h. Te 3-(4,5-dimethylthiazol-2-yi)-2,5diphenyltetrazolium bromide (MTT) was used to assess cell proliferation, and the absorbance was detected by Tecan A-5082 sunrise (Austria) at the same time points from day 1 to day 5. Te colony information assay was performed by seeding 5 * 10 2 cells in 6-well plates after transfection for 48 h. After being cultured for 2-3 weeks, the colonies were fxed and stained.
2.11. Transwell Assay. Te transwell assay was performed by seeding 5 × 10 4 cells in the upper chambers with 10% FBS medium, while the lower chambers were added with 20% FBS medium. After incubating for 10-14 h., the cells in the upper chamber were fxed and stained by a three-step set (Termo Scientifc, USA). Images were taken by a light microscope (Olympus, Japan) at 100 * magnifcation, the migrated cells were counted and calculated.
2.12. Western Blot. RIPA bufer (Solarbio, China) with PMSF (Termo Scientifc, USA) was utilized to extract cellular proteins, and the concentration was detected using a BCA kit (Termo Scientifc, USA). Te 10% SDS-PAGE gels were used for protein separation and then the separated proteins were transferred to PVDF membranes. Te 5% skimmed milk block was applied for one hour and then incubated at 4°C overnight with diluted primary antibodies. After three times washing with TBST, the membranes were incubated with secondary antibodies for one hour at room temperature. Finally, the blots were detected by ECL (Millipore, USA). Te antibodies information was listed in the supplementary fle: Table S2. 2.13. Statistical Analysis. Mean and standard deviation (SD) were utilized for presenting data. Te signifcance of differences between the experimental and control groups was determined by Student's t-test. Results with P < 0.05 were considered with statistical signifcance. Te calculations of all data were conducted using IBM SPSS software (version 22.0, USA).

Identifcation and Functional Analysis of DEGs.
To identify diferentially expressed genes related to TME, the specimens were divided into high-and low-score groups according to stromal score, immune score, and ESTIMATE score. Finally, 987 DEGs were identifed according to stromal score, containing 914 upregulated genes and 73 down-regulated genes (Figures 1(a), 1(c) and 1(d)). Correspondingly, 1267 DEGs were obtained according to the immune score, including 954 upregulated genes and 313 downregulated genes (Figures 1(b)-1(d)). Te Venn plot was used to identify 788 upregulated genes and 65 downregulated genes in both the immune and stromal components (Figures 1(c) and 1(d)). Ten, these 853 DEGs were considered TME-related DEGs. To further investigate the biological functions and capable pathways associated with these TME-related DEGs, the analyses of KEGG and GO were performed. Te results of KEGG analysis shows that DEGs were involved in some immune-related activities, for instance, cytokine-cytokine receptor interaction and chemokine signaling pathway (Figure 2(a)). While the results of GO analysis, which contains molecular functions (MF), biological processes (BP), and cellular components (CC), indicated that the TME-related DEGs were primarily enriched in immune-related functions, including GO: 0042110 (T-cell activation), GO:0001772 (immunological synapse), and GO:0140375 (immune receptor activity)

PPI Network and Univariate Cox Regression Analysis.
After obtaining TME-related DEGs through the above analyses, the PPI network had been constructed using the String dataset to investigate the interactions among these genes and the result was visualized using Cytoscape ( Figure 3(a)). Ten, the univariate Cox regression analysis was performed to explore the correlation between DEGs and the prognosis of patients with TC. Four genes were identifed as risk factors of TC prognosis (Figure 3(b)), which were set as the main characters of our following analysis. To further investigate their functions, GSEA was conducted among these four DEGs, which indicated that they were largely engaged in immune-related events, such as the JAK-STAT signaling pathway, Tcell receptor pathway, and autoimmune thyroid diseases. Besides, it is worth mentioning that specifc enrichment results also showed the signifcant role of these genes in cancer, including the P53 signaling pathway, pathways in cancer, and TGF-β signaling pathway (Figures 3(c)-3(f )).

Te Role of 4
Core DEGs in TME of TC. Ten, the infltrating data of 22 diferent types of immune cells in tumor tissues was identifed using the CIBERSORT algorithm (Figures 4(a) and 4(b)). Te diferent types of TIICs were discovered to be closely involved in varying core gene expression in the TME of TC cells, and detailed results were established as violin plots (Figure 4(c)). Te interrelation between the proportions of TIICs and core DEGs expression was shown as the scatter plots, respectively, including LRRN4CL ( Figure 5 , which showed these 4 core genes were positively correlated with dendritic cells (DCs), monocytes, (MC) and neutrophils, while negatively associated with NK cells, plasma cells, CD8+ T cells, and M0 macrophages. Terefore, these results indicated that LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8 play a crucial role in the immune activities of TC cells.

Te Correlation between 4 Core DEGs and Clinical
Signatures. In the previous study, the core DEGs were identifed as LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8. Te TC specimens were divided into high-and lowgroups according to the expression of core DEGs, then we investigated the role of these core genes in the clinical outcome of TC patients, which showed that all these 4 DEGs were related to poor OS in patients with TC (P � 0.042, P � 0.008, P � 0.036, and P � 0.046, respectively) ( Figure 6(d)). Additionally, the expressions of these DEGs were relatively associated with the N stage, especially for CAPN8 (P � 3.1e − 05) ( Figure 6(c)). Te Wilcoxon rank-sum test indicated the expression levels of 4 DEGs among normal and tumor tissue in paired or unpaired samples, which indicated HS3ST3A1 and CAPN8 expressed signifcantly higher in tumor tissues, while the other two genes, LRRN4CL and PCOLCE2, were not (Figures 6(a) and 6(b)). Terefore, HS3ST3A1 and CAPN8 were selected for further investigation to verify their malignancy behavior in experiments.

HS3ST3A1 and CAPN8 Were Associated with Higher
Tumor Stage. Te GEPIA (https://gepia.cancer-pku.cn) was employed for further investigation of the correlation between HS3ST3A1 or CAPN8 with clinical stages, which results showed that both HS3ST3A1 and CAPN8 signifcantly differed among stages (Figure 7(a)). In our previous study, the expression of DEGs seemed to be inextricably linked to the N stage. Terefore, we collected 10 paired samples of thyroid cancer with lymph node metastasis (LNM) and nonlymph node metastasis (No LNM), respectively. RT-qPCR results indicated that the mRNA levels of HS3ST3A1 and CAPN8 were remarkably higher in tumor tissues compared with normal tissues. Furthermore, their mRNA expression levels in tumors with LNM were notably higher compared with those without LNM (Figure 7(b)). Ten, we collected 30 cases of TC patient tissues, including normal, tumors with LNM and tumors with no LNM. Te IHC results indicated that HS3ST3A1 and CAPN8 expression levels were both higher in tumors with LNM than in tumors with no LNM or normal tissues (Figure 7(c)). Tese fndings indicated that HS3ST3A1 and CAPN8 were intimately correlated to higher clinical stages and lymph node metastasis.

Downregulation of HS3ST3A1 or CAPN8 Inhibit Malignancy Behaviors of Tyroid Cancer Cells.
To verify the cellular function of HS3ST3A1 and CAPN8, we transfected shHS3ST3A1 and shCAPN8 into papillary thyroid cancer cells BCPAP and TPC-1, which downregulation was detected by Western blot (Figure 7(d)). Furthermore, cells with       shHS3ST3A1 and shCAPN8 showed less proliferation rates (Figure 7(e)) and generated fewer colonies (Figure 7(f )) than the control groups. As for the transwell assay, cells with shHS3ST3A1 and shCAPN8 also behaved less aggressively (Figure 7(g)). Tese results demonstrated that the downregulation of HS3ST3A1 and CAPN8 was able to inhibit the proliferation and invasion of thyroid cancer cells in vitro.

Discussion
It has been reported in a huge number of studies that TME plays critical roles in multiple cancer types and promotes the progression of molecular classifcation systems and treatment strategies [15]. Our study identifed LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8 as TME-related genes associated with prognosis in TC samples collected from the TCGA database, and verifed two of them, HS3ST3A1 and CAPN8, through a clinical specimen and cellular functional experiments. LRRN4CL lacks in-depth research, while only one study reported that the upregulation of the cell surface protein LRRN4CL promoted metastases of pulmonary in mice [16]. HS3ST3A1 encodes the enzyme 3-O-sulfotransferase, which catalyzes the biosynthesis of 3-O-sulfated heparan sulfate, a specifc subtype of heparan sulfate (HS). It has been reported that HS3ST3A1 is involved in respiratory papillomatosis [17] and human immunodefciency virus (HIV) infection [18]. Besides, HS3ST3A1 was found to be upregulated in lung cancer tissue compared with normal lung tissue and associated with the progression of lung cancer [19]. PCOLCE2 is reported to encode a functional collagenbinding protein procollagen C-proteinase enhancer (PCPE2) [20], and evidence has proved that PCOLCE2 was able to perform as a biomarker for prognostic prediction in colorectal cancer patients [21]. CAPN8 belongs to the calpain family and exhibits restricted expression patterns [22] and has been proposed to be involved in vesicle trafcking [23]. However, there are no investigations on the role of these four genes in thyroid cancer. Synthesizing previous studies, there are hardly any reports on their functions in the tumor microenvironment, but their roles as prognosis indicators are certainly clarifed in several types of cancer. It is our study that collects LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8 all together for the frst time and sets them as biomarkers related to TME in thyroid cancer, which might become potential candidate targets for TC immunotherapy.
It is commonly recognized that infammation and autoimmunity are risk factors for TC [24]. Evidence also indicated that TC patients might beneft from targeting cancer-related infammation, which provided new strategies for diagnosis and treatment [25]. Studies of immune infltration in TC have made many major advances, particularly in the study of primary immune cells, such as NK cells, TAMs, MCs, DCs, CD8+ T cells, B cells, neutrophils,  Journal of Oncology and Tregs. [24,26,27]. In our study, the proportion of TIICs was estimated using the CIBERSORT algorithm, which indicated that TME-related DEGs had a notable association with specifc immune cells, including monocytes, neutrophils, T cells, and so on. Previous studies revealed that TAMs, MCs, DCs, Tregs, and neutrophils were positively related to TC progression [28][29][30][31][32]. While NK cells, iDCs, CD8+ T cells, and B cells are negatively associated with TC     progression [33][34][35][36][37][38]. However, some studies clarifed that immune cells are related to various outcomes in cancer [39]. It is sure that further studies will be conducted to explore the role of diferent immune cells plays in cancer for better prognostic evaluation.
As for TME-related genes in thyroid cancer, a previous study has identifed 30 hub genes by constructing the PPI network and set CXCL10 as the top hub gene [14]. Our study has diferent logical methods and conducted further data mining. First, we obtained more gene expression and clinical profles so that our study might have more credibility. CXCL10 was a specifc diferential expressed TME-related gene in the PPI network with the most nodes, but it tended to lose its signifcance when performing univariate Cox regression analysis. Furthermore, we analyzed the immune infltrating profles of TC and concluded the correlation between immune cells and DEGs using the CIBERSORT algorithm. Tese make sense for better identifying TMErelated biomarkers and prognostic indicators in patients with TC.
HS3ST3A1 and CAPN8 were considered to serve as prognostic predictors in our study. Clinical specimens were collected and subjected to RT-qPCR and IHC staining, which showed higher HS3ST3A1 and CAPN8 expression levels in tumor tissues, especially in tumors with LNM. In vitro, the downregulation of HS3ST3A1 and CAPN8 presented the inhibition of proliferation and invasion in papillary thyroid cancer cells. Tese fndings indicated that TME-related genes HS3ST3A1 and CAPN8 function in the immune process and contribute to tumor development.
However, it is necessary to acknowledge the limitations of this study. Te biases were unavoidable because the data were primarily obtained from the TCGA database. Besides, we chose two of the genes for experimental verifcation; the other two genes, LRRN4CL and PCOLCE2 were excluded due to their adverse expression both in paired and unpaired tissues. As for clinical specimens, our data lacks survival information due to the good outcome of TC patients. Furthermore, our experiments mainly focused on papillary thyroid cancer cells, which is the most common subtype of TC, but the verifcation in other thyroid cancer subtypes is lacking. Animal experiments and the depth of molecular mechanisms still need further investigation.
Overall, we used the ESTIMATE algorithm to identify DEGs related to the TME in TC specimens obtained from the TCGA dataset. LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8 were observed as potential prognostic indicators for patients with TC. Furthermore, HS3ST3A1 and CAPN8 were highly expressed in thyroid tumor tissue, especially in tumors with LNM, and their downregulation can inhibit the proliferation and invasion of thyroid cancer cells. However, the underlying molecular mechanisms of tumor micrometastasis and the potential clinical value for early diagnosis still require further experimental study.

Conclusion
In summary, we identifed several TME-related DEGs in TC, among them, LRRN4CL, HS3ST3A1, PCOLCE2, and CAPN8 were remarkably involved in the regulation of the immune activities in the TME and poor clinical outcomes. Additionally, the malignancy behaviors of HS3ST3A1 and CAPN8 in the tumor process were verifed through tissue and cellular experiments. Our fndings indicated that HS3ST3A1 and CAPN8 served as prognostic biomarkers of TC and might bring new insights into the development of efective therapeutic strategies for patients with TC.

Data Availability
All data generated or analyzed during this study are included in this article and its supplementary fles.

Consent
Tis study was approved by the Institutional Review Board of Tianjin Medical University Cancer Institute and Hospital. All methods were performed in accordance with the relevant guidelines and regulations. Voluntary written informed consent had been obtained from each participant.