ORIGINAL RESEARCH article

Front. Mol. Biosci., 08 March 2021

Sec. Molecular Diagnostics and Therapeutics

Volume 7 - 2020 | https://doi.org/10.3389/fmolb.2020.602068

Identification of HOXA1 as a Novel Biomarker in Prognosis of Head and Neck Squamous Cell Carcinoma

  • Department of Otolaryngology, The First Affiliated Hospital of Bengbu Medical College, Bengbu, China

Abstract

Hox genes, a highly conserved homolog in most animals, play vital functions in cell development and organ formation. In recent years, researchers have discovered that it can act as a tumor regulator, and its members can participate in tumorigenesis by regulating receptor signaling, cell differentiation, apoptosis, migration, EMT, and angiogenesis. Hox genes and which major members play a vital role in the progress of head and neck squamous cell carcinoma (HNSCC) is still unclear. After analyzing the expression differences and prognostic value of all Hox genes through the TCGA-HNSC database, we use histochemistry stains in 52 pairs of HNSCC slices to verify the expression level of the key member-HOXA1. In correlation analysis, we found that high HOXA1 expression is related to poor pathological grade (p = 0.0077), advanced T stage (p = 0.021) and perineural invasion (PNI) (p = 0.0019). Furthermore, we used Cox univariate and multivariate regression analysis to confirm the independent predictive power of HOXA1 expression. To explore the underlying mechanisms behind HOXA1, we ran GSVA and GSEA and found fourteen mutual signaling pathways, including neuroprotein secretion and transport, tumor-associated signaling pathways, cell adhere junction and metabolic reprogramming. Finally, we found that the high expression of HOXA1 is significantly related to the decrease of CD8+ T cell infiltration and the decline of DNA methylation level. Our findings demonstrated that HOXA1, as a notable member of the HOX family, maybe an independent prognostic indicator in HNSCC.

Introduction

Head and neck cancer was one of the most common malignant upper respiratory tract tumors globally, of which squamous cell carcinoma (HNSCC) accounts for about 90% (Bray et al., 2018). Although the continuous improvement of treatment such as surgery, radiotherapy, chemotherapy, targeted therapy and immunotherapy, the 5-year survival rate of HNSCC is still stagnant (Cancer Genome Atlas Network, 2015). Senile patients are mainly those who are habitually ingesting alcohol and tobacco, while HPV infection has led to an increase in young patients recently, especially oropharyngeal squamous cell carcinoma (Chaturvedi et al., 2011). Among them, patients with HPV-negative suffer from irreversible damage due to long-term smoking and drinking, who are more tolerant of the treatments with a worse prognosis.

Hox genes is a subgroup of the homeobox transcription factor family, which aggregates in 4 types of chromosomal loci (A-D), with 39 genes in humans. Hox gens plays its main function by encoding transcription factors or directly targeting downstream genes as monomers or homodimers (Alexander et al., 2009). Members of Hox genes can participate in cell development, organ formation and even carcinogenesis by regulating receptor signaling, cell differentiation, apoptosis, migration, epithelial-mesenchymal transition (EMT) and angiogenesis (Cho et al., 2012). Hox genes is tissue-specific and may play a key role in developmental regulation in one tissue type (Economides and Capecchi, 2003), while in another tissue, it is an accomplice of tumor malignant biological behavior (Wang et al., 2007). Previous studies have found that many Hox genes are abnormally expressed in hematological malignancies and various solid tumors. There are inadequate studies on the involvement of Hox genes in the occurrence and development of HNSCC. Therefore, studying the role of Hox genes in HNSCC is of great significance.

In this study, we used the method of bioinformatics mining database of TCGA-HNSC combined with immunohistochemical (IHC) verification to find a novel prognostic biomarker, HOXA1 (Homeobox A1) in Hox genes. HOXA1 is the first Hox gene to be discovered that can be upregulated during carcinogenesis by gastrin (Hofsli et al., 2005). In breast cancer, HOXA1 has been identified as a breast epithelial oncogene whose forced expression is sufficient to induce the transformation of immortalized human mammalian epithelial cells into aggressive cancer cells (Zhang et al., 2006; Mohankumar et al., 2008). In previous reports, the abnormal expression of HOXA1 is related to the prognosis of breast cancer, lung cancer, liver cancer, prostate cancer, oral cancer, gastric cancer and melanoma (Shah and Sukumar, 2010). We evaluated the abnormal expression of HOXA1 in HNSCC and confirmed its independent predictive power. At the same time, we analyzed the correlation between the expression level and clinicopathological characteristics. According to the expression level, we ran gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) to further analyze the biological pathways of HNSCC pathogenesis related to the expression of HOXA1. Finally, results of TIMER, CIBERSORT and MEXPRESS analysis indicated that HOXA1 was correlated with CD8+ T cell infiltration and DNA methylation level.

Materials and Methods

HOX Genes Expression Datasets and Relevant Clinicopathological Features

HNSCC dataset (528 cases) contain RNA-Seq (HTSeq) and clinicopathological characteristics from the cancer genome atlas (TCGA-HNSC) were downloaded (https://cancergenome.nih.gov/). Screened out 499 cases, including the expression of 39 Hox genes and corresponding clinicopathological features are used in the subsequent analysis. R Package “survminer” were used to calculate the cutoff value of Hox genes. The 499 cases were divided into a high-expression group (n = 396) and a low-expression group (n = 103) according to the cutoff value of HOXA1 expression (overall survival, OS: 1.0525, disease-free survival, DFS: 2.9749). Clinical features of HNSCC patients had been shown in Table 1.

Table 1

Clinical characteristicsTotal
(n = 499)
Expression of HOXA1p-Value*
HighLow
(n = 396)(n = 103)
Age, n (%)0.8741
    <65311 (62.3)248 (62.6)63 (61.2)
    ≥65188 (37.7)148 (37.4)40 (38.8)
Gender, n (%)0.0086
    Female133 (26.7)110 (27.8)23 (22.3)
    Male366 (73.3)286 (72.2)80 (77.7)
Histologic grade, n (%)0.0023
    G1-2362 (72.5)300 (75.8)62 (60.2)
    G3-4121 (24.2)84 (21.2)37 (35.9)
    NA16 (3.2)12 (3.0)4 (3.9)
AJCC Stage, n (%)1.0000
    I-II114 (22.8)91 (23.0)23 (22.3)
    III-IV371 (74.3)296 (74.7)75 (72.8)
    NA14 (2.8)9 (2.3)5 (4.9)
T stage, n (%)0.0080
    T1-2176 (35.3)129 (32.6)47 (45.6)
    T3-4308 (61.7)258 (65.2)50 (48.5)
    NA15 (3.0)9 (2.3)6 (5.8)
N classification, n (%)0.1761
    N0 (negative)238 (47.7)197 (49.7)41 (39.8)
    N1-3(positive)239 (47.9)185 (46.7)54 (52.4)
    NA22 (4.4)14 (3.5)8 (7.8)
M classification, n (%)1.0000
    M0469 (37.3)374 (94.4)95 (92.2)
    M15 (0.2)4 (1.0)1 (1.0)
    NA25 (50.3)18 (4.5)7 (6.8)
HPV_STATUS_P16, n (%)1.0000
    Positive34 (6.8)22 (5.6)12 (11.7)
    Negative78 (15.6)61 (15.4)17 (16.5)
    NA387 (75.56)313 (79.0)74 (71.8)
ANGIOLYMPHATIC_ INVASION, n (%)0.5906
    YES218 (43.7)92 (23.2)28 (27.2)
    NO120 (24.0)174 (43.9)44 (42.7)
    NA161 (32.3)130 (32.8)31 (30.1)
PERINEURAL_INVASION, n (%)0.00197
    YES185 (37.1)142 (35.9)23 (22.3)
    NO165 (33.1)133 (33.6)52 (50.5)
    NA149 (29.8)121 (30.6)28 (27.2)
Smoking pack / years, n (%)0.5314
    ≥40161 (32.3)127 (32.1)34 (33.0)
    <40126 (25.3)104 (26.3)22 (21.4)
    NA212 (42.4)165 (41.7)47 (45.6)
ALCOHOL_HISTORY_ DOCUMENTED, n (%)1.0000
    YES331 (66.3)262 (66.2)69 (67.0)
    NO157 (31.5)125 (31.5)32 (31.1)
    NA11 (2.2)9 (2.3)2 (1.9)
Primary Tumor Site, n (%)0.4934
    Oropharynx378 (75.8)296 (74.7)82 (79.6)
    Larynx111 (22.2)91 (23.0)20 (19.4)
    Hypopharynx10 (2.0)9 (2.3)1 (1.0)

Clinical features of TCGA-HNSC patients (n = 499).

*

The chi-square test. Bold values indicate statistically significant, P < 0.05.

Differentially Expressed Genes and Clinical Pathological Parameters Analysis

First, the R package “limma” was used to analyze the differential expression of 39 Hox genes between tumor and normal samples in 499 patients and differentially expressed genes (DEGs) with |log2 fold change (FC)| > 1 and adjusted P-values < 0.05 were selected for subsequent analysis. The prognostic genes with the predictive ability of OS and DFS were screened out by the analysis of the Kaplan–Meier survival curve and the log-rank test (P < 0.05). The correlation analysis between the expression of HOXA1 and clinicopathological characteristics were evaluated by the R package “combn.” Furthermore, Cox univariate and multivariate regression analysis were performed on the individual clinical variables [age (years ≥65/ <65), gender (male/female), grade (G1&2/G3&4), AJCC Stage (I&II/III&IV), T stage (T1–2/T3–4), N classification (N0/N1–3), M classification (M0/M1), HPV_STATUS_P16 (positive/negative), angiolymphatic invasion (Yes/No), perineural invasion (Yes/No), smoking (packs/year ≥40/ <40), alcohol history documented (Yes/No) and primary tumor site (Oropharynx/ Larynx/Hypopharynx)]with the HOXA1 expression value and high or low-expression level. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated using the R package “coxph.” The chi-square test was used to assess the statistical significance of expression value between the two groups.

Analysis of Immune Infiltration and DNA Methylation

The fraction of twenty-two immune cell types of each TCGA-HNSC sample was yielded through R package “CIBERSORT” by estimating relevant subsets of RNA profiles (https://cibersort.stanford.edu/) (Supplementary Material 1). Only samples with P < 0.05 were conceived to perform the further analysis of comparing differential the tumor-infiltrating immune cells (TICs) level between the high and low-HOXA1 level group. In addition, another immune infiltration data is derived from TIMER (Tumor Immune Estimation Resource, https://cistrome.shinyapps.io/timer/) and DNA methylation data came from MEXPRESS database (https://mexpress.be/).

Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA)

The R package “GSVA” was used to evaluate the signaling pathways (5529 gene sets) and Gene Ontology (GO) terms (10,192 gene sets) in the TCGA-HNSC database (Hänzelmann et al., 2013). Then, the R package “limma” was used to analyze the significant differential gene sets between the two groups. Besides, 499 HNSCC cases were divided according to the median expression value of HOXA1. GSEA was then performed to detect the gene sets that were enriched in the gene rank for identifying potential the pathways most associated with HOXA1. The annotated gene sets of c2.cp.kegg. v7.1. symbols in the Molecular Signatures Database (MSigDB) was selected in GSEA version 4.1.0 (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). We performed 1,000 times of permutations and filtered the results by normalized enrichment scores (NES), nominal p-value and false discovery rate (FDR) q-value (Subramanian et al., 2005).

Tissue Specimen Collection and Immunohistochemistry (IHC)

The study was approved by the Clinical Research Ethics Committee of The First Affiliated Hospital of Bengbu Medical College. Written informed consent was obtained from all-volunteer patients. In all, 52 pairs of HNSCC paraffin sections were collected from The First Affiliated Hospital of Bengbu Medical College (Supplementary Material 2). In the progress of antigen retrieval, the sections were placed in 10 mM citrate buffer solution (pH = 6.0) and boiling in a pressure cooker for 15 min. After cooling to 37°C, 3% H2O2 was added and incubated for 10 min to deplete the activity of endogenous peroxidase. After blocking by Goat serum, the sections were stained with an anti-HOXA1 polyclonal antibody (1:100, Thermo Fisher Scientific, catalog # PA5-36164, RRID AB_2553401) for a night at 4°C-refrigerator. After soaking and rinsing thoroughly with PBS three times, the sections were stained with secondary antibodies. Finally, the positive expression of the HOXA1 protein was observed by a microscope after adding substrate and hematoxylin staining. The expression density was scored by two experienced pathologists, following the single-blind trial technique. HOXA1 protein staining was assessed by the sum of the percentage of positive nuclear cells: − (0–25%), + (26–50%), ++ (51–100%). At last, we evaluated HOXA1 expression into two levels, low (− to +) and high (++). Fifty-two pairs of IHC pictures were taken for each slide, and images (200×) were analyzed using software ImageJ (US National Institute of Health, Bethesda, Maryland, USA, https://imagej.nih.gov/).

Statistical Analysis

Statistical analyses were conducted with software R (v3.6.0) and GraphPad Prism 7 (GraphPad Software, California, USA). Differences among variables were analyzed with two-tailed Student's t-test and chi-square test. All significance tests were carried out at the 0.05 level.

Results

Differentially Expression and Prognostic Value of HOX Genes in HNSCC

The analysis flow chart of this article is shown in Figure 1. Most Hox genes (30/39, Supplementary Material 3) have significant differences in expression between 499 tumor samples and 43 tumor-adjacent normal samples (Figure 2A). Among DEGs, only HOXA1 have predictive ability both in OS and DFS (Supplementary Figure 1), and worse clinical outcomes in patients with high expression of HOXA1 (Figures 2B–E).

Figure 1

Figure 2

HOXA1 Is Upregulated in HNSCC

In TCGA-HNSC, HOXA1 expression in tumor tissues was significantly higher than adjacent tissues (p = 1.241e-10) (Figure 3A), the differential analysis of HOXA1 expression in 43 paired samples got the same result (p = 9.984e−09) (Figure 3B). We performed IHC in 52 pairs of HNSCC paraffin sections for further verification (Figures 3E–G). HOXA1 protein immunostaining was mainly found in the cell nucleus but widely distributed in cancer nest. Among them, HOXA1 protein was detected in 45 cases (86.5%, + to ++) of tumor samples, while HOXA1 protein was detected in only ten sections (19.2%, + to ++) of adjacent pieces. There were 38 cases with low positive staining (− to +) and 14 patients with high positive staining (++) among the tumor samples. Our results also showed that HNSCC samples have a higher mean percentage of HOXA1-positive area (Figure 3C) and HOXA1-positive cells (Figure 3D) than paracancerous samples (P < 0.0001). The numeric dimension of the areas (mean ± SD, median) used for the IHC evaluation were shown in Table 2.

Figure 3

Table 2

The numeric dimensionAdjacent tissues
(mean ± SD, median)
HNSCC tissues
(mean ± SD, median)
P-value
HOXA1-positive area2.86 ± 1.29, 2.9410.44 ± 3.24, 10.48<0.0001
HOXA1-positive cells16.69 ± 7.60, 15.5040.87 ± 11.62, 42.00<0.0001

The numeric dimension of the areas (mean ± SD, median) used for the IHC evaluation.

Association of HOXA1 Expression With Clinicopathological Features

The expression of HOXA1 was highly statistically significant correlated with pathology grade (G1&2 vs. G3&4, p = 0.0077), local invasion (T1–2 vs. 3–4, p = 0.023), perineural invasion (No vs. Yes, p = 0.0019) and HPV infection (positive vs. negative, p = 0.001) (Figures 4A–L). The expression of HOXA1 in oropharyngeal cancer was significantly lower than that in laryngeal cancer (p = 0.021), and there was no significant difference between the other two tumor locations (Figure 5A). We conducted further survival analysis in three tumor locations and found that the expression of HOXA1 only has prognostic value in oropharyngeal cancer (P < 0.05), while no significant predictive value in hypopharyngeal and laryngeal cancer (p > 0.05) in TCGA-HNSC patients (Figures 5B–D).

Figure 4

Figure 5

Independent Predictive Power of HOXA1 Expression Value and Expression Level

In Cox univariate regression analysis, the HRs of high-expression vs. low-expression of HOXA1 and expression value of HOXA1 for overall survival (OS) was 2.316 (P < 0.001, CI: 1.515–3.540) and 1.236 (P < 0.001, CI: 1.115–1.369), respectively (Figure 6A). In further multivariate Cox regression analysis, the expression level and expression value of HOXA1 still have a good predictive ability for OS [expression level (HR = 1.927, p = 0.010, CI: 1.173–3.166) and expression value (HR = 1.144, p = 0.047, CI: 1.002–1.307)] (Figure 6B). The results of the above analysis indicated HOXA1 was an independent prognostic factor, and the expression level of HOXA1 seem better.

Figure 6

GSVA and GSEA Reveal HOXA1-Related Pathways and Molecular Functions

To explore the potential carcinogenic function of HOXA1, we used both GSVA and GSEA both approaches. After scoring 499 patients by the GSVA, the significant differences in signaling pathways and GO functions associated with the expression level of HOXA1 were obtained, as shown in Figures 7A,B. Besides, we ran GSEA to figure out the significant differential signaling pathways related to the expression value of HOXA1 (Table 3), including neuroprotein secretion and transport, the tumor-associated signaling pathway, cell adhere junction and metabolic reprogramming (Figures 8A–D). After the intersection of the results of GSVA and GSEA, 14 mutual pathways mainly focusing on the following four aspects: neurotrophin, vesicular transport, protein export; Ubiquitin mediated proteolysis, Glycan biosynthesis; gap junction, adherens junction, focal adhesion, regulation of actin cytoskeleton; Pathways in cancer including WNT/ERBB/ NOTCH/TGF-beta.

Figure 7

Table 3

MSigDB collectionGene set nameNESNOM p-valFDR q-val
c2.cp.kegg.v7.1.symbols.gmtKEGG_WNT_SIGNALING_PATHWAY1.96100.075
KEGG_GLYCOSYLPHOSPHATIDYLINOSITOL_GPI_ANCHOR_BIOSYNTHESIS1.9520.0020.054
KEGG_ERBB_SIGNALING_PATHWAY1.7600.0030.137
KEGG_UBIQUITIN_MEDIATED_PROTEOLYSIS1.8290.0040.100
KEGG_N_GLYCAN_BIOSYNTHESIS2.0570.0040.038
KEGG_ADHERENS_JUNCTION1.8750.0050.079
KEGG_PURINE_METABOLISM1.6790.0060.141
KEGG_OTHER_GLYCAN_DEGRADATION1.7580.0080.125
KEGG_REGULATION_OF_ACTIN_CYTOSKELETON1.7120.0150.135
KEGG_PROTEIN_EXPORT1.7200.0190.138
KEGG_NOTCH_SIGNALING_PATHWAY1.6640.0260.132
KEGG_GAP_JUNCTION1.6140.0280.155
KEGG_TGF_BETA_SIGNALING_PATHWAY1.6740.0290.137
KEGG_FOCAL_ADHESION1.6960.0310.132
KEGG_INOSITOL_PHOSPHATE_METABOLISM1.6290.0340.161
KEGG_PATHWAYS_IN_CANCER1.5910.0380.146
KEGG_NEUROTROPHIN_SIGNALING_PATHWAY1.6230.0390.160
KEGG_SNARE_INTERACTIONS_IN_VESICULAR_TRANSPORT1.5880.0450.142

Gene sets enriched related with expression of HOXA1.

Gene sets with NOM P-value < 0.05 and FDR q-value < 0.25 are considered as significant.

MsigDB, Molecular Signatures Database; NES, normalized enrichment score; NOM, nominal; FDR, false discovery rate.

Figure 8

Immune Infiltration and DNA Methylation Correlating With HOXA1 Expression

More and more studies have found that some oncogenes significantly affect the tumor immune microenvironment, including infiltration of B cell, CD8+ T cell, CD4+ T cell, Macrophage, Neutrophil, Dendritic cell, etc. The HOXA1 expression in the TIMER database is associated with a decrease of B cell / CD8+T cell infiltration and an increase of CD4+ T cell infiltration in total HNSCC patients, a reduction of B cell and CD8+T cell infiltration in HPV(+) HNSCC patients, and a decrease of CD8+T cell infiltration and an increase of CD4+ T cell infiltration in HPV(−) HNSCC patients (Figures 9A–C). From the results of immune infiltration in total and HPV subtype HNSCC patients, it is found that the high expression of HOXA1 significantly leads to the reduction of CD8+ T cell infiltration. Besides, we performed the CIBERSORT analysis to screen the differential TICs of the high and low-HOXA1 expression groups and found that the results were similar to the TIMER analysis results. The CIBERSORT analysis results showed that the proportion of M0 macrophages was positively correlated with HOXA1 expression level, and the ratio of CD8+T cells, naïve B cells, CD4 memory activated T cells and follicular helper T cells were negatively associated with the level of HOXA1 expression (Figure 9D).

Figure 9

The correlation analysis between DNA methylation level and HOXA1 expression showed that high HOXA1 expression with less promoter methylation in tumor samples (Pearson correlation coefficients ranging from −0.166 to −0.528 for promoter region probes) in the MEXPRESS database (Figure 10). We also found that the probes such as cg03116258, cg07450037 and cd12686016 of the seventh chromosome had the strongest correlation in DNA methylation of HOXA1, which was worthy of further exploration and research.

Figure 10

Discussion

There are only sporadic reports in the literature on the role of HOX genes in HNSCC. Jinyun Li et al. used bioinformatics to analyze the HOXA genes' role in laryngeal squamous cell carcinoma (LSCC), in which HOXA10, HOXA11 and HOXA13 expression and clinical characteristics (M stage and gender) are independent prognostic biomarkers for the prognosis of LSCC patients (Li et al., 2020). The expression of HOXA members in LSCC is negatively correlated with the level of methylation. HOXA7, HOXA10 and HOXA13, which are highly expressed in LSCC tissues, have a weak positive correlation with tumor purity. However, there is no significant correlation between the TICs and HOXA members. Studies in oral squamous cell carcinoma (OSCC) showed that HOXA genes are more enriched than normal mucosal tissues. High expression of HOXA1 can promote OSCC proliferation and maybe an independent prognostic indicator for patients (Bitu et al., 2012). Both AMOT-p80 and miR-146a, identified as HOXD10 targets in HNSCC, can represent therapeutic targets for specific tumor stages (Hakami et al., 2014). High expression of HOXD10 promotes tumor proliferation and migration, while low expression of HOXD10 can promote tumor invasion and metastasis by regulating AMOT-p80. In this paper, we showed the expression profile of 39 members of the HOX genes in HNSCC, as well as the most important independent prognostic indicator and potential biological role of HOXA1 in HNSCC, providing novel insights for further research on HOXA1 as a potential target of HNSCC.

The expression of HOXA1 is tissue-specific (Fagerberg et al., 2014) (https://www.proteinatlas.org), and as tumor modulators, it may cover a wide range of functionality in types of tumor dependent on the cellular context (Shah and Sukumar, 2010). HOXA1 played a role of an oncogene in breast carcinoma (Liu et al., 2019), lung carcinoma (Abe et al., 2006), oral squamous cell carcinoma (Bitu et al., 2012), hepatocellular carcinoma (Zha et al., 2012), melanoma (Wardwell-Ozgo et al., 2014), gastric cancer (Yuan et al., 2016), leukemia (Chen et al., 2019) and retinoblastoma (RB) (Lyv et al., 2020). Inversely, it is shown as anti-oncogenes in pancreatic carcinoma (Ohuchida et al., 2012) and small cell lung cancer (SCLC) (Xiao et al., 2014). Nilay et al. pointed out that the abnormal expression of Hox genes triggers oncogenesis or tumor suppression, which may be caused by her mechanism of action (Shah and Sukumar, 2010). On the level of expression, Hox genes are associated with oncogenesis when expressed at an increased level in tumor tissues not normally seen in that tissue type. In the results of our HNSCC IHC, the expression of HOXA1, which is mainly located in the cell nucleus, is extremely low in adjacent normal tissues while it is widely distributed in cancer nest. Thus, it can be seen that the abnormal up-regulation of HOXA1 in HNSCC, especially in the nucleus, indicates that HOXA1 mainly plays a role in oncogenesis.

E-cadherin activated Rac signaling pathway through the transcription target of HOXA1 (Zhang et al., 2006), which can promote breast cancer cell survival by upregulating cyclinD1, c-Myc, and BCL-2 (Mohankumar et al., 2007, 2008). HOXA1 also takes effect independent of transcriptional regulation (Taminiau et al., 2016). In proteomics research, it is found that HOXA1 also participates in cell signal transduction as a functional protein, such as pathways of focal adhesion, receptor tyrosine kinase function, TGF/BMP and TNF/NF-kB (Lambert et al., 2012). In breast cancer, HOXA1 transforms immortalized human breast epithelial cells stimulated by p44/42 MAPK pathway, signal transducer and activator of transcription 3 (STAT3) and STAT5B into an invasive malignant phenotype. WNT, TGF-beta and FGF pathways may synergistically affect HOXA1 expression, and WNT7A directly maintains expression of HOXA1, making for lung cancer recurrence (Calvo et al., 2000). HOXA1 dedifferentiates melanoma cells into a highly aggressive cell state accompanied by TGF activation (Wardwell-Ozgo et al., 2014). Down-regulation of HOXA1 can increase the expression of E-cadherin while reducing the expression of Snail and MMP-3, which can inhibit ERK1/2 and AKT to suppress the growth, migration, invasion and metastasis of prostate cancer (Wang et al., 2015).

The methylation of DNA is inversely proportional to transcription activity, which also affects the genome's stability. Yeon et al. evaluated the methylation status of 15 genes of promoter CpG island, including HOXA1, that are involved in breast cancer progression (Park et al., 2012). The expression level of HOXA1 in subtypes of breast cancers are related to the disordered epigenetic activity, in which the methylation frequency of Luminal B is the highest (< 0.001). Meanwhile, HOXA1 hypermethylation is related to the absence of HER2 neu expression and proliferation of breast cancer cells (P < 0.05) (Pilato et al., 2013). Through the MEXPRESS database, we found that more HOXA1 expression in HNSCC samples with less promoter methylation, which possibly foreshadowing the lower DNA methylation level of oncogenes, may lead to the occurrence of tumors.

In addition to DNA methylation, many non-coding RNAs (ncRNAs) as another way of apparent epigenesis have been found to participate in the control of gene expression and activation of pathways, which in turn participates in the development of cancer and other diseases. Many studies have shown that circRNAs, lncRNAs and miRNAs play a key role in tumorigenesis, invasion, metastasis and drug resistance through the function of competitive RNA (ceRNA). Here we summarize Non-coding genes act as regulators of HOXA1 in cancers (Table 4).

Table 4

Non-coding formNon-coding RNAsCancer typesFunctionsRef. (PMID)
microRNAmiR-193a-5pBreast cancerTumor promotor32497022
miR-100Lung cancerChemotherapy resistance24559685
Nasopharyngeal CarcinomaGrowth and proliferation32021301
miR-218Lung cancerGefitinib Resistance32084702
Cervical cancerProliferation, colony formation, migration and invasion30896864
miR-181bSCLCRegulatory of p53 signaling pathway29658571
miR-181a-5pMMApoptosis resistance29228867
miR-216b-5pCervical cancerProliferation31114990
miR-210SCLC/melanomaImmune escape22962263
miR-145Oral squamous cell carcinomaCancer inhibition31138758
miR-10aPCTumor promotor22407312
miR-30cMalignant giant-cell tumor of boneMetastasis and growth29164581
miR-30bEsophageal cancerGrowth, migration and invasion28189678
miR-99aHCCMigration and invasion31186723
miR-99 familyEpithelial cellProliferation and migration24312487
miR-515RBTumor promotor32561925
LncRNALncRNA SNHG1Breast cancerTumor promotor32497022
LncRNA CCAT1Lung cancerGefitinib Resistance32084702
MMApoptosis resistance29228867
LINC00152Cervical cancerProliferation31114990
lncRNA HOTAIRSCLCMultidrug resistance26707824
LncRNA HOTAIRM1Breast cancerAvoiding silence of HOXA132284737
Endometrial CancerProliferation, Migration, Invasion and EMT31853186
GBMGrowth and invasion30376874
LncRNA ZFPM2-AS1RBTumor promotor32561925
CircRNACircEIF4G2Cervical cancerProliferation, colony formation, migration and invasion30896864

Non-coding genes act as regulators of HOXA1 in cancers.

SCLC, small cell lung cancer; RB, retinoblastoma; GBM, glioblastoma multiforme; PC, Pancreatic cancer; MM, Multiple Myeloma; HCC, hepatocellular carcinoma cells; EMT, epithelial–mesenchymal transition.

Tumor immunity research is in full swing, and we are still full of unknowns on how the abnormal expression of HOXA1 affects immune function. In the study of endometrial cancer, high expression of HOXA1 in Myeloid-derived suppressor cells (MDSCs) promotes tumor progression and hinders anti-tumor immune response (Li et al., 2019). Tumors adapt to hypoxic stress and create an immunosuppressive microenvironment suitable for survival (Keith et al., 2011). Studies have found that in hypoxic areas of SCLC and melanoma, miR-210 regulates the expression of HOXA1 and inhibits the anti-tumor effect of cytotoxic T lymphocyte (Noman et al., 2012). Regarding the results of TIMER and CIBERSORT analysis, we suspected that the high expression of HOXA1 in HNSCC might reduce CD8+ T cell infiltration, leading to an immunosuppressive environment not conducive to patient's survival outcomes.

Previous studies of pancreatic cancer and glioma showed that the close interaction between tumors and neurons is not conducive to the clinical prognosis. The density of newborn nerve cells in the tumor microenvironment is often associated with tumor growth and highly invasive biological behaviors (Magnon et al., 2013; Venkataramani et al., 2019). The latest study by Moran Amit et al. found that the loss of TP53 in HNSCC can promote tumor growth by affecting neuronal remodeling (Amit et al., 2020). So far, we have found that the function of p53 in SCLC and osteosarcoma is related to the regulation of HOXA1 (Zhao et al., 2015; Zhang et al., 2018). In mouse development studies, it was found that at embryonic stage (D7-9), HOXA1 expression precisely regulates the development of neural crest precursor cells, which are essential for tumor microenvironment remodeling (Gavalas et al., 2001). In this paper, we found that high HOXA1 expression is related to PNI by mining TCGA-HNSC data. According to the analysis of GSVA and GESA, we speculated that HNSCC might act on neurons by secreting exosomes containing neurotrophic proteins. However, this result lacked evidence concerning that no cases that meet PNI were found after careful screening of the paraffin blocks by pathologists. Whether HNSCC hijacks surrounding neurons to promote tumor growth through the regulation of HOXs is still worth exploring.

In summary, HOXA1 may mainly exert tumor regulation functions through the following four mechanisms: spatiotemporal expression disorder (high expression in cancer tissues, especially in the cell nucleus), abnormal protein-signaling pathways (neuronal protein secretion), mediation of immunosuppression (inhibition of CD8+ T cell infiltration) and aberrant epigenetic activity (high expression of HOXA1 leads to a decrease in DNA methylation level). Our study demonstrated that HOXA1 is highly expressed in HNSCC, which may be an independent prognostic factor in TCGA-HNSC patients, especially for oropharyngeal cancer patients. At present, HOXA1 small interfering RNA (siRNA) nanoparticles have achieved initial success in mouse models of breast cancer and are effective (75%) to reduce the incidence of tumors (Brock et al., 2014). More random cohort studies and mechanism experiments will help confirm the regulatory effects of HOXA1 on HNSCC.

Statements

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/repository.

Ethics statement

The studies involving human participants were reviewed and approved by Clinical Research Ethics Committee of The First Affiliated Hospital of Bengbu Medical College. The patients/participants provided their written informed consent to participate in this study.

Author contributions

The project was proposed by HL and SM. HL provided paraffin tissue sections and wrote the manuscript. MW assisted in editing the manuscript. JZ downloaded and mined the relevant data of TCGA. XW and MZ did immunohistochemistry and assisted in organizing data. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by the Anhui Provincial Key Research and Development Project (202004j07020007), Anhui Provincial Natural Science Foundation (1608085MH236), Key Project of Natural Science of Bengbu Medical College (BYKY2019086ZD), and Science Research Project of Bengbu Medical College (BYKF1855).

Acknowledgments

Thanks to Zhaoyi Lu (Department of Otolaryngology-Head and Neck Surgery, Xiangya Hospital, Central South University) for suggestions on the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2020.602068/full#supplementary-material

Supplementary Figure 1

Kaplan-Meier survival analysis of (A) overall survival and (B) disease free survival) according to the expression of all HOX genes in TCGA-HNSC patients. The log-rank test was used to calculate the p-value.

Supplementary Material 1

CIBERSORT results of the fraction of twenty-two immune cell types of each TCGA-HNSC sample.

Supplementary Material 2

Clinicopathological characteristics of 52 HNSCC patients collected from The First Affiliated Hospital of Bengbu Medical College.

Supplementary Material 3

Results of differential expression analysis in all HOX genes.

References

  • 1

    AbeM.HamadaJ.-I.TakahashiO.TakahashiY.TadaM.MiyamotoM.et al. (2006). Disordered expression of HOX genes in human non-small cell lung cancer. Oncol. Rep.15, 797802. 10.3892/or.15.4.797

  • 2

    AlexanderT.NolteC.KrumlaufR. (2009). Hox genes and segmentation of the hindbrain and axial skeleton. Annu. Rev. Cell Dev. Biol.25, 431456. 10.1146/annurev.cellbio.042308.113423

  • 3

    AmitM.TakahashiH.DragomirM. P.LindemannA.Gleber-NettoF. O.PickeringC. R.et al. (2020). Loss of p53 drives neuron reprogramming in head and neck cancer. Nature578, 449454. 10.1038/s41586-020-1996-3

  • 4

    BituC. C.Destro MFdS. S.CarreraM.da SilvaS. D.GranerE.KowalskiL. P.et al. (2012). HOXA1 is overexpressed in oral squamous cell carcinomas and its expression is correlated with poor prognosis. BMC Cancer12:146. 10.1186/1471-2407-12-146

  • 5

    BrayF.FerlayJ.SoerjomataramI.SiegelR. L.TorreL. A.JemalA. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin.68, 394424. 10.3322/caac.21492

  • 6

    BrockA.KrauseS.LiH.KowalskiM.GoldbergM. S.CollinsJ. J.et al. (2014). Silencing HoxA1 by intraductal injection of siRNA lipidoid nanoparticles prevents mammary tumor progression in mice. Sci. Transl. Med.6:217ra2. 10.1126/scitranslmed.3007048

  • 7

    CalvoR.WestJ.FranklinW.EricksonP.BemisL.LiE.et al. (2000). Altered HOX and WNT7A expression in human lung cancer. Proc. Natl. Acad. Sci. U.S.A.97, 1277612781. 10.1073/pnas.97.23.12776

  • 8

    Cancer Genome Atlas Network (2015). Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature517, 576582. 10.1038/nature14129

  • 9

    ChaturvediA. K.EngelsE. A.PfeifferR. M.HernandezB. Y.XiaoW.KimE.et al. (2011). Human papillomavirus and rising oropharyngeal cancer incidence in the United States. J. Clin. Oncol.29, 42944301. 10.1200/JCO.2011.36.4596

  • 10

    ChenS.-L.QinZ.-Y.HuF.WangY.DaiY.-J.LiangY. (2019). The role of the gene family in acute myeloid leukemia. Genes10:621. 10.3390/genes10080621

  • 11

    ChoH.-S.ToyokawaG.DaigoY.HayamiS.MasudaK.IkawaN.et al. (2012). The JmjC domain-containing histone demethylase KDM3A is a positive regulator of the G1/S transition in cancer cells via transcriptional regulation of the HOXA1 gene. Int. J. Cancer131, E179E189. 10.1002/ijc.26501

  • 12

    EconomidesK. D.CapecchiM. R. (2003). Hoxb13 is required for normal differentiation and secretory function of the ventral prostate. Development130, 20612069. 10.1242/dev.00432

  • 13

    FagerbergL.HallströmB. M.OksvoldP.KampfC.DjureinovicD.OdebergJ.et al. (2014). Analysis of the human tissue-specific expression by genome-wide integration of transcriptomics and antibody-based proteomics. Mol. Cell. Proteomics MCP13, 397406. 10.1074/mcp.M113.035600

  • 14

    GavalasA.TrainorP.Ariza-McNaughtonL.KrumlaufR. (2001). Synergy between Hoxa1 and Hoxb1: the relationship between arch patterning and the generation of cranial neural crest. Development128, 30173027.

  • 15

    HakamiF.DardaL.StaffordP.WollP.LambertD. W.HunterK. D. (2014). The roles of HOXD10 in the development and progression of head and neck squamous cell carcinoma (HNSCC). Br. J. Cancer111, 807816. 10.1038/bjc.2014.372

  • 16

    HänzelmannS.CasteloR.GuinneyJ. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics14:7. 10.1186/1471-2105-14-7

  • 17

    HofsliE.ThommesenL.YadetieF.LangaasM.KusnierczykW.FalkmerU.et al. (2005). Identification of novel growth factor-responsive genes in neuroendocrine gastrointestinal tumour cells. Br. J. Cancer92, 15061516. 10.1038/sj.bjc.6602535

  • 18

    KeithB.JohnsonR. S.SimonM. C. (2011). HIF1α and HIF2α: sibling rivalry in hypoxic tumour growth and progression. Nat. Rev. Cancer12, 922. 10.1038/nrc3183

  • 19

    LambertB.VandeputteJ.RemacleS.BergiersI.SimonisN.TwizereJ.-C.et al. (2012). Protein interactions of the transcription factor Hoxa1. BMC Dev. Biol.12:29. 10.1186/1471-213X-12-29

  • 20

    LiJ.YeM.ZhouC. (2020). Expression profile and prognostic values of family members in laryngeal squamous cell cancer. Front. Oncol.10:368. 10.3389/fonc.2020.00368

  • 21

    LiX.PangL.YangZ.LiuJ.LiW.WangD. (2019). LncRNA HOTAIRM1/HOXA1 axis promotes cell proliferation, migration and invasion in endometrial cancer. Onco. Targets. Ther.12, 1099711015. 10.2147/OTT.S222334

  • 22

    LiuJ.LiuJ.LuX. (2019). HOXA1 upregulation is associated with poor prognosis and tumor progression in breast cancer. Exp. Ther. Med.17, 18961902. 10.3892/etm.2018.7145

  • 23

    LyvX.WuF.ZhangH.LuJ.WangL.MaY. (2020). Long noncoding RNA ZFPM2-AS1 knockdown restrains the development of retinoblastoma by modulating the microRNA-515/HOXA1/Wnt/β-catenin axis. Invest. Ophthalmol. Vis. Sci.61:41. 10.1167/iovs.61.6.41

  • 24

    MagnonC.HallS. J.LinJ.XueX.GerberL.FreedlandS. J.et al. (2013). Autonomic nerve development contributes to prostate cancer progression. Science341:1236361. 10.1126/science.1236361

  • 25

    MohankumarK. M.PerryJ. K.KannanN.KohnoK.GluckmanP. D.EmeraldB. S.et al. (2008). Transcriptional activation of signal transducer and activator of transcription (STAT) 3 and STAT5B partially mediate homeobox A1-stimulated oncogenic transformation of the immortalized human mammary epithelial cell. Endocrinology149, 22192229. 10.1210/en.2007-1320

  • 26

    MohankumarK. M.XuX. Q.ZhuT.KannanN.MillerL. D.LiuE. T.et al. (2007). HOXA1-stimulated oncogenicity is mediated by selective upregulation of components of the p44/42 MAP kinase pathway in human mammary carcinoma cells. Oncogene26, 39984008. 10.1038/sj.onc.1210180

  • 27

    NomanM. Z.BuartS.RomeroP.KetariS.JanjiB.MariB.et al. (2012). Hypoxia-inducible miR-210 regulates the susceptibility of tumor cells to lysis by cytotoxic T cells. Cancer Res.72, 46294641. 10.1158/0008-5472.CAN-12-1383

  • 28

    OhuchidaK.MizumotoK.LinC.YamaguchiH.OhtsukaT.SatoN.et al. (2012). MicroRNA-10a is overexpressed in human pancreatic cancer and involved in its invasiveness partially via suppression of the HOXA1 gene. Ann. Surg. Oncol.19, 23942402. 10.1245/s10434-012-2252-3

  • 29

    ParkS. Y.KwonH. J.ChoiY.LeeH. E.KimS.-W.KimJ. H.et al. (2012). Distinct patterns of promoter CpG island methylation of breast cancer subtypes are associated with stem cell phenotypes. Modern Pathol.25, 185196. 10.1038/modpathol.2011.160

  • 30

    PilatoB.PintoR.De SummaS.LamboR.ParadisoA.TommasiS. (2013). HOX gene methylation status analysis in patients with hereditary breast cancer. J. Hum. Genet.58, 5153. 10.1038/jhg.2012.118

  • 31

    ShahN.SukumarS. (2010). The Hox genes and their roles in oncogenesis. Nat. Rev. Cancer10, 361371. 10.1038/nrc2826

  • 32

    SubramanianA.TamayoP.MoothaV. K.MukherjeeS.EbertB. L.GilletteM. A.et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A.102, 1554515550. 10.1073/pnas.0506580102

  • 33

    TaminiauA.DraimeA.TysJ.LambertB.VandeputteJ.NguyenN.et al. (2016). HOXA1 binds RBCK1/HOIL-1 and TRAF2 and modulates the TNF/NF-κB pathway in a transcription-independent manner. Nucleic Acids Res.44, 73317349. 10.1093/nar/gkw606

  • 34

    VenkataramaniV.TanevD. I.StrahleC.Studier-FischerA.FankhauserL.KesslerT.et al. (2019). Glutamatergic synaptic input to glioma cells drives brain tumour progression. Nature573, 532538. 10.1038/s41586-019-1564-x

  • 35

    WangH.LiuG.ShenD.YeH.HuangJ.JiaoL.et al. (2015). HOXA1 enhances the cell proliferation, invasion and metastasis of prostate cancer cells. Oncol. Rep.34, 12031210. 10.3892/or.2015.4085

  • 36

    WangZ.DahiyaS.ProvencherH.MuirB.CarneyE.CoserK.et al. (2007). The prognostic biomarkers HOXB13, IL17BR, and CHDH are regulated by estrogen in breast cancer. Clin. Cancer Res.13, 63276334. 10.1158/1078-0432.CCR-07-0310

  • 37

    Wardwell-OzgoJ.DogrulukT.GiffordA.ZhangY.HeffernanT. P.van DoornR.et al. (2014). HOXA1 drives melanoma tumor growth and metastasis and elicits an invasion gene expression signature that prognosticates clinical outcome. Oncogene33, 10171026. 10.1038/onc.2013.30

  • 38

    XiaoF.BaiY.ChenZ.LiY.LuoL.HuangJ.et al. (2014). Downregulation of HOXA1 gene affects small cell lung cancer cell survival and chemoresistance under the regulation of miR-100. Eur. J. Cancer50, 15411554. 10.1016/j.ejca.2014.01.024

  • 39

    YuanC.ZhuX.HanY.SongC.LiuC.LuS.et al. (2016). Elevated HOXA1 expression correlates with accelerated tumor cell proliferation and poor prognosis in gastric cancer partly via cyclin D1. J. Exp. Clin. Cancer Res.35:15. 10.1186/s13046-016-0294-2

  • 40

    ZhaT.-Z.HuB.-S.YuH.-F.TanY.-F.ZhangY.ZhangK. (2012). Overexpression of HOXA1 correlates with poor prognosis in patients with hepatocellular carcinoma. Tumour Biol.33, 21252134. 10.1007/s13277-012-0472-6

  • 41

    ZhangX.EmeraldB. S.MukhinaS.MohankumarK. M.KraemerA.YapA. S.et al. (2006). HOXA1 is required for E-cadherin-dependent anchorage-independent survival of human mammary carcinoma cells. J. Biol. Chem.281, 64716481. 10.1074/jbc.M512666200

  • 42

    ZhangY.LiX.-J.HeR.-Q.WangX.ZhangT.-T.QinY.et al. (2018). Upregulation of HOXA1 promotes tumorigenesis and development of non-small cell lung cancer: a comprehensive investigation based on reverse transcription-quantitative polymerase chain reaction and bioinformatics analysis. Int. J. Oncol.53, 7386. 10.3892/ijo.2018.4372

  • 43

    ZhaoS.NiuF.XuC.-Y.YeL.BiG.-B.ChenL.et al. (2015). Microarray and ChIP-seq data analysis revealed changes in p53-mediated transcriptional regulation in Nutlin-3-treated U2OS cells. Mol. Med. Rep.12, 42844290. 10.3892/mmr.2015.3933

Summary

Keywords

hox genes, homeobox A1, head and neck squamous cell carcinoma, perineural invasion, prognosis, biomarker, tumor

Citation

Li H, Wang X, Zhang M, Wang M, Zhang J and Ma S (2021) Identification of HOXA1 as a Novel Biomarker in Prognosis of Head and Neck Squamous Cell Carcinoma. Front. Mol. Biosci. 7:602068. doi: 10.3389/fmolb.2020.602068

Received

02 September 2020

Accepted

12 November 2020

Published

08 March 2021

Volume

7 - 2020

Edited by

Satyendra Chandra Tripathi, All India Institute of Medical Sciences Nagpur, India

Reviewed by

José Manuel Lopes, Universidade Do Porto, Portugal; Ajay Dixit, University of Minnesota Twin Cities, United States; Bin Zhao, Xiamen University, China

Updates

Copyright

*Correspondence: Shiyin Ma

This article was submitted to Molecular Diagnostics and Therapeutics, a section of the journal Frontiers in Molecular Biosciences

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics