HIF1A overexpression predicts the high lymph node metastasis risk and indicates a poor prognosis in papillary thyroid cancer

Objective To investigate the value of Hypoxia-inducible factor 1 A (HIF1A) in predicting lymph node metastasis (LNM) stage and clinical outcomes of papillary thyroid cancer (PTC) patients. Materials and methods The HIF1A gene expression analysis in PTC was performed by bioinformatics approaches followed by evaluating its protein level using immunohistochemistry analysis. The role of HIF1A in predicting the LNM stage was evaluated by logistic regression analysis, nomogram construction, and receiver operating characteristic (ROC) analysis. We performed survival analyses to determine its prognostic value. Enrichment analysis was conducted, and immune cell infiltration and stromal content were evaluated to examine the underlying mechanism of HIF1A in PTC. Results HIF1A transcription and protein levels were significantly high in PTC tissue (P < 0.05). Its overexpression predicted high LNM risk and unfavorable prognosis for PTC patients (P < 0.05). Cox regression analysis revealed HIF1A as an independent prognostic biomarker for the disease-free interval (DFI) (P < 0.01). In addition, HIF1A was positively related to tumor-suppressive immunity but was negatively correlated with anti-tumor immunity. HIF1A upregulation was also associated with increased stromal content. Conclusions HIF1A overexpression is an independent predictor for worse DFI in PTC. The HIF1A expression may affect the prognosis of PTC patients through immune- and stroma-related pathways. Our study provides new insight into the role of HIF1A in PTC biology and clinical management.


Introduction
Thyroid cancer is the most common endocrine malignancy with an increasing incidence globally in recent years. Thyroid cancer ranked fourth in the number of all estimated new cases of cancers in 2021, while it was ranked fifth and sixth in 2020 and 2019, respectively [1][2][3]. It includes papillary, follicular, medullary, and undifferentiated types [4]. Papillary thyroid cancer (PTC) is the most prevalent histological type of carcinoma with an 80% prevalence [5]. PTC patients in early-stage have an overall 5-year survival of over 90% using established methods [6]. Although most PTC patients have a favorable prognosis, a small proportion of patients still suffer from local invasion or distant metastasis, leading to a low survival time [7]. Lymph node metastasis (LNM) occurs in 14%-64% of PTC patients and is associated with recurrence [8]. Besides, the enhanced angiogenesis pathway also results in the progression of PTC [9]. Therefore, it is urgent to find novel biomarkers for predicting the LNM stage and survival to improve the prognosis of PTC patients.
Hypoxia-inducible factor 1 (HIF-1) contains an oxygen-regulated subunit (HIF1A) and a constitutively expressed subunit [10]. Upon hypoxia, HIF-1 binds to the hypoxia response elements (HREs), regulating multiple pathological and physiological processes including apoptosis, glucose metabolism, and angiogenesis [11]. Specifically, HIF1A participates in adapting cancer cells to lower oxygen levels and it may be a late recurrence biomarker [12]. The activation of HIF1A promotes the Warburg effect by switching from oxidative phosphorylation to glycolysis [11]. HIF1A expression is low in most cells under normal conditions, but HIF1A is often significantly elevated under hypoxia [13]. HIF1A is overexpressed in cancers with involvement in biological processes of tumor cell survival, angiogenesis, metastasis, and tumor therapy [14,15]. Besides, HIF1A is expressed in lymphocytes, dendritic cells, neutrophils, and macrophages, thereby regulating the immune microenvironment [16]. The previous study has demonstrated upregulated HIF1A expression in breast tumors compared with normal breasts and served as a potential prognostic biomarker for cancer recurrence [12]. However, there is no systemic research on the expression, clinical significance, and underlying mechanism of HIF1A in PTC.
This study aimed to find a potent biomarker to predict LNM stage and prognosis in PTC. In the area of big data, bioinformatics as a powerful tool for data mining was adopted for HIF1A expression analysis. We also examined value of HIF1A in predicting LNM stage and prognosis of PTC patients. The enrichment analysis, evaluation of immune cell infiltration, and stromal content were employed to elucidate the underlying mechanism of HIF1A in PTC.

Data acquisition
Due to the small number of normal thyroid samples in the TCGA database, the data from TCGA combined with GTEx (including 504 thyroid cancer samples and 338 normal thyroid samples) were obtained. The transcriptome profiles (HTSeq-FPKM), phenotype information, and survival data were downloaded from the GDC TCGA sets in the UCSC Xena database (https://xenabrowser.net/). Inclusion criteria: PTC samples; patients had complete survival and expression information. At last, we included 497 PTC samples for following analyses, and they were divided into low-and high-HIF1A expression groups according to the median HIF1A expression. The patient characteristics were shown in Table S1. Besides, as the following inclusion criteria: >100 samples; containing thyroid cancer and normal samples; with expression data, the GSE33630 dataset (containing 59 tumorous and 46 normal samples) was generated to verify HIF1A expression.

Expression analysis of HIF1A
TIMER (http://timer.comp-genomics.org/) was used to explore the differential HIF1A mRNA level in various cancers. Then, HIF1A gene expression in PTC was analyzed using the GDC TCGA together with GTEx data, validated by the GSE33630 dataset. Following the manufacturer's protocol, we carried out immunohistochemistry (IHC) analysis for acquiring its protein expression.
Next, the relationship between HIF1A levels and clinical parameters of PTC patients including age, gender, LNM stage, and subtype was evaluated using the chi-square test, and univariate logistic regression analysis. Female was set as a reference level for gender, N0 stage for the LNM stage, and PTC-classical for the histological subtype.

The clinical significance of HIF1A
Firstly, the association of the LNM stage with HIF1A expression and the other clinical parameters was determined using logistic regression analysis. After that, we constructed a nomogram. Besides, we drew the ROC curves to examine the value of the HIF1A together with clinical factors in predicting the LNM stage.
Following this, the effect of HIF1A expression on survival was analyzed using the Kaplan-Meier plotter method. The survival differences between two expression groups were compared by the log-rank test. Further, the independent predictors for PTC were identified using Cox regression analysis.

Functional enrichment analysis
The limma package was employed for differential expression analysis between the two HIF1A groups. The significant differentially expressed genes (DEGs) were identified by setting |log 2 foldchange |>1 and P < 0.05 as the selection criterion. Gene ontology annotations and possible pathways enriched by the significant DEGs were analyzed with clusterProfiler package [17].

Gene set enrichment analysis (GSEA)
To further explore the biological signaling pathways that HIF1A was involved in, GSEA was executed between two HIF1A expression groups. The gene set was permutated 1000 times and HIF1A expression level was used as a phenotypic label.

Single-sample gene-set enrichment analysis (ssGSEA)
Moreover, we defined the enrichment level of pathways in PTC samples as the ssGSEA score [18]. The gene set represents the collection of all marker genes of a pathway. higher HIF1A mRNA level in PTC than that in normal group using the TCGA combined with GTEx data, validated by the GSE33630 dataset. (C) The higher protein level of HIF1A in PTC tissue using immunohistochemistry analysis. *P < 0.05, ***P < 0.001.

Evaluation of the immune cell infiltration and stromal content
CIBERSORT algorithm is a tool for estimating the abundances of 22 immune cells using transcriptomic data [19,20]. We used this algorithm to calculate relative immune cell proportions in two HIF1A expression groups. In addition, ESTIMATE algorithm was used to calculate the Immune Score (represents the infiltration of immune cells in tumor tissue), the Stromal Score (captures the presence of stroma in tumor tissue), and the ESTIMATE Score (equaled the sum of Stromal Score and Immune Score) [21].

Statistical analysis
SPSS software (version 23.0) and R software (version 4.0.3) were used for all statistical analyses. Differences in each two-group comparison and among at least three groups were determined by a t-test and one-way ANOVA, respectively. Pearson's correlation test was adopted for correlation analysis. P < 0.05 was considered statistically significant.

HIF1A expression is increased in PTC patients
The HIF1A mRNA level in human cancers was visualized in Fig. 1A. Then, we found that HIF1A gene was elevated in PTC compared to normal thyroid (P < 0.001). GSE33630 data further confirmed the HIF1A overexpression in the tumor tissue (P < 0.001) (Fig. 1B). IHC staining exhibited that HIF1A protein was obviously upregulated in PTC (Fig. 1C). These results indicated high HIF1A expression in PTC.

HIF1A is linked to LNM stage and subtype
We analyzed the distribution of different clinicopathological parameters in two HIF1A groups. There was no significantly different distribution of age and gender (P > 0.05), while the N1 stage was mainly distributed in the high HIF1A group (P < 0.001). Besides, the histological subtype was significantly related to HIF1A expression (P < 0.001) ( Table 1).
When considering HIF1A as a continuous variable, HIF1A was not associated with age and gender as well ( Fig. 2A-B). However, HIF1A was significantly linked to LNM stage and subtype (P < 0.001) ( Fig. 2C-D). As expected, patients at the N1 stage had higher HIF1A expression compared with those at the N0 stage (P < 0.05) (Fig. 2C). Additionally, PTC-tall cell had the highest HIF1A expression, followed by PTC-classical, and PTC-follicular (Fig. 2D).
Through the logistic regression analysis, age and gender had no remarkable correlation with the HIF1A expression. The high LNM stage was associated with HIF1A overexpression with an odds ratio (OR) of 2.071 (P < 0.001). Besides, PTC-tall cell and PTC-follicular were also significantly correlated with HIF1A expression (P < 0.05) ( Table 2).

HIF1A overexpression predicted the high LNM risk
Due to the significant relation between HIF1A expression and the LNM stage, and the prognostic value of the LNM stage in PTC [22], we further analyzed its value in predicting LNM stage. HIF1A, age, gender, and subtype could independently predict the LNM stage (P < 0.05). Of note, HIF1A overexpression was connected with high LNM stage risk with an OR of 1.796 (P < 0.001) ( Table 3). After integrating independent factors into nomogram construction, HIF1A contributed the most to the LNM risk (Fig. 3A). Moreover, the AUC of the ROC curve of HIF1A was 0.701, which was superior to age (AUC = 0.571), gender (AUC = 0.534), and subtype (AUC =   Abbreviations: LNM stage, lymph node metastasis; PTC, papillary thyroid cancer.

Table 3
The association of the LNM stage with the HIF1A expression and clinicopathological parameters using logistic regression analysis. 0.611), indicating that HIF1A presented an acceptable performance in predicting the LNM stage (Fig. 3B).

High HIF1A mRNA expression leads to poor prognosis
We plotted the survival curves to explore the prognostic significance of HIF1A in PTC. Fig. 4A-D shows that patients with high HIF1A expression tended to have a shorter DFI and PFI, but had no significant difference in DSS and OS. The stratified analysis presented that HIF1A had a significant impact on DFI and PFI in PTC patients aged below 55 years, females, and the PTC-classical subtype (P < 0.05). In addition, HIF1A had a close relation with DFI in those at N1 status but significantly affected PFI in those at N0 status (P < 0.05) ( Table 4).
To identify independent risk factors for DFI and PFI, Cox regression analysis was conducted. In the univariate analysis, the LNM stage and HIF1A were significantly related to DFI and PFI (P < 0.05). Next, these two significant factors were included in the multivariate analysis. They could independently predict DFI but not PFI (Table 5).

HIF1A participates in immune-and stromal-associated pathways
To elucidate biological function of HIF1A for PTC, we first screened 1062 downregulated and 764 upregulated significant DEGs between two expression groups. Then, these DEGs were used for functional enrichment analysis to explore the HIF1A-related pathways. The top 21 significant terms of BP, CC, and MF enrichment analysis are presented in Fig. 5A-C. Notably, in terms of BP, HIF1A was mainly enriched in the regulation of the immune system process, biological adhesion, regulation of immune response, adaptive immune response, humoral immune response, B cell-mediated immunity, and complement activation. Additionally, the major KEGG pathways in which they were involved were cytokine-cytokine receptor interaction, cell adhesion molecules, hematopoietic cell lineage, focal adhesion, ECM-receptor interaction, and autoimmune thyroid disease (Fig. 5D).
To further explore the molecular mechanisms affected by HIF1A in PTC, GSEA was executed using GDC TCGA data. JAK-STAT signaling pathway, cell adhesion molecules, NOD-like receptor signaling pathway, focal adhesion, chemokine signaling pathway, Toll-like receptor signaling pathway, and ECM-receptor interaction were enriched in the high HIF1A expression group (Fig. 6). The combined results uncovered that hematopoietic cell lineage, JAK-STAT signaling pathway, Toll-like receptor signaling pathway, cell adhesion molecules, focal adhesion, and ECM-receptor interaction might be important pathways involved in PTC.
Following this, the enrichment levels of hematopoietic cell lineage, JAK-STAT signaling pathway, Toll-like receptor signaling pathway, cell adhesion molecules, focal adhesion, and ECM-receptor interaction were analyzed using the ssGSEA. The HIF1A gene had a significantly positive relationship with the enrichment levels of the six pathways (all P < 0.001) (Fig. 7A-F). Therefore, HIF1A might affect the PTC progression mainly through positive regulation of immune-related pathways including hematopoietic cell lineage, JAK-STAT, and Toll-like receptor, and by stromal-associated pathways including cell adhesion molecules, focal adhesion, and ECM-receptor interaction.

Association of HIF1A with immune cell infiltration and stromal content
We have demonstrated that HIF1A was involved in pathways linked to the immune and stroma; hence, we evaluated the immune cell infiltration and stromal content. The abundance of most immune cells was significantly different between high-and low-HIF1A expression groups (Fig. 8A). HIF1A was negatively correlated with T cells CD8, NK cells activated, monocytes, and macrophages M2, but positively associated with dendritic cells activated, and T cells regulatory (Fig. 8B). These results indicated that HIF1A expression exhibited negative correlations with anti-tumor immunity, and presented a positive association with tumor-suppressive immunity.
Subsequently, ESTIMATE algorithm was used for evaluation of stromal content. HIF1A was positively related to the immune score, stromal score, and ESTIMATE score (Fig. 9A-C); indicating that the elevated HIF1A level was linked to a more abundant stromal content in cancer, another potential mechanism by which HIF1A regulated PTC progression.

Discussion
PTC is the most common subtype of thyroid cancer, and patients with PTC always exhibit a good prognosis with radical surgery and radiation therapy; nevertheless, some patients still have aggressive disease and can develop distant metastasis [23]. ANXA10 has been proven to be upregulated in PTC tissue, and its high expression independently predicted the poor prognosis for PTC patients [24].  CXCL10 was a potential core gene of the PTC microenvironment and was significantly associated with prognosis and immune cells in PTC [25]. Despite the rapid progress in PTC biomarkers, the OS rate of PTC patients has not significantly improved in the past 10 years [26]. Our study demonstrates the role of HIF1A in PTC and reveals its associated regulatory pathways by performing a comprehensive analysis of open-access databases.
Hypoxia is a common condition found in a range of solid tumors, which plays an essential role in tumor occurrence and development [27]. Tumor cells adapt to a hypoxic environment mainly regulated by HIF-1 [28]. As a subunit of HIF-A, HIF1A has been implicated as a tumor promoter gene, which was overexpressed in tumors of the ovarian, bladder, and prostate [29][30][31]. Moreover, HIF1A expression was highly expressed in thyroid cancer compared with that in normal thyroid [32]. However, there is a lack of systemic research about the expression, clinical significance, and underlying mechanism of HIF1A in PTC. Our study also showed an upregulation of HIF1A in PTC tissues compared with the normal thyroid tissue, which was validated by IHC analysis. Besides, HIF1A Abbreviations: HR, hazard ratio; 95% CI, 95% confidence interval; LNM, lymph node metastasis; PTC, papillary thyroid cancer. expression was significantly higher in PTC patients at the N1 stage, and the HIF1A overexpression was significantly related to high LNM risk with an OR = 1.796, which was similar to the result of Liu's study [33]. Moreover, we found that HIF1A was superior to other clinical factors in predicting the LNM stage. It has been demonstrated that HIF1A upregulates epithelial-mesenchymal transition-related transcription factors [34]. These results suggested that HIF1A might represent a crucial role in the invasion, metastasis, and progression of PTC. The prognostic value of HIF1A in PTC was also assessed using GDC TCGA data. Patients with high HIF1A expression had shorter DFI and PFI but had no significant difference in DSS or OS compared to those with low HIF1A expression. Further, Cox regression analysis  revealed HIF1A as a prognostic factor for poor DFI. We have demonstrated that patients with high HIF1A expression are more likely to develop lymph node metastases, which may explain the unfavorable prognosis in the high HIF1A expression group. HIF1A binds to the hypoxia response elements in the promoter region of target genes and is mainly involved in adaptive changes, enabling tumor cells to survive and proliferate in a hypoxic environment, thereby promoting malignant phenotypes and aggressive tumor behaviors [35]. Therefore, the authors hypothesized that HIF1A overexpression might contribute to poor clinical outcomes by providing a growth advantage under hypoxic stress, accelerating tumor invasion and metastasis.
Subsequently, we explored the mechanism of HIF1A regulating PTC. HIF1A was mainly involved in immune-associated pathways of JAK-STAT, Toll-like receptor, hematopoietic cell lineage, as well as stroma-associated pathways including cell adhesion molecules, focal adhesion, and ECM-receptor interaction pathway. Aberrant activation of the Toll-like receptor signaling pathway leads to NF-κB signaling activation and overexpression of inflammatory cytokines such as IL-6, IL-1β and TNF-α [36]. JAK-STAT signaling mediates almost all immune regulatory processes, including those involved in tumor cell recognition and tumor-driven immune escape [37]. The presence of stimulatory or inhibitory signals governs the innate and adaptive immune activity, controlling effective immune surveillance and facilitating escape. These signals consist of various antigens, foreign ligands, nucleic acids, as well as cytokines and secreted factors that are induced by IFN, TGF-β, and NF-κB pathways regulated by the activation of JAK-STAT [38]. The ECM is a highly dynamic and complex molecular network that surrounds cells in tissues. Focal adhesion is a subcellular structure that provides strong adhesion to the ECM and serves as a scaffold for many integrin-involved signaling pathways [39]. In addition, focal adhesion represents an essential role in tumor cell invasion and cell migration [40]. CIBERSORT analysis also showed that HIF1A was negatively correlated with T cells CD8, NK cells activated, monocytes, and macrophages M2, but positively associated with dendritic cells activated, and T cells regulatory. High expression of T cells CD8 indicates obvious T cell reaction in thyroid cancer patients and its infiltration contributed to better OS. Dendritic cells were increased in thyroid cancer tissue, suppressing the immune response [41]. Macrophages M2 favor the Th2 immune response, promoting tumor progression [42]. Taken together, HIF1A expression exhibited negative correlations with anti-tumor immunity and presented a positive association with tumor-suppressive immunity. In addition to tumor cells and immune cells, stromal cells are also important components in the tumor microenvironment. Cytokines secreted by stromal cells activate signaling pathways regulating immune infiltration and EMT [43]. ECM components provide biochemical and biomechanical cues for cells, which are vital in breast cancer progression and metastasis [44]. The focal adhesion pathway has also been revealed to participate in the EMT process, and its activation can alter cell glycolysis, and induce drug resistance [45,46]. Thus, HIF1A upregulation might lead to a worse prognosis through the activation of these oncogenic pathways.
The comprehensive study is potentially valuable in advancing the understanding of the relationship between HIF1A and PTC, but some limitations still exist. Although we performed the IHC analysis to determine the HIF1A protein expression in tissues of PTC and normal thyroid, the role of HIF1A in tumor metastasis and immune cell infiltration require further study in the future. Second, most of the analyses performed in the study were based on HIF1A mRNA level, and deeper analysis based on protein level would increase the evidence level of the data. These results should be validated in a large-scale population in the future.
In conclusion, the HIF1A expression was increased in PTC tissue, and its expression was positively linked to the LNM stage. HIF1A overexpression could predict high LNM risk and poor DFI, indicating that HIF1A might serve as a reliable biomarker for PTC patients. Moreover, HIF1A might influence survival of PTC patients via regulating immune-as well as stroma-related pathways.

Disclosure of interest
The authors report no conflict of interest.

Consent for publication
Not applicable.

Data sharing statement
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board of the First Hospital of Lanzhou University.