Differentially Expression of LHXs in HNSCC Tumor Samples
The flow diagram of this study is shown in Figure 1. We first detected 72 differentially expressed genes (DEGs) in the homeobox gene superfamily between HNSCC tumor and normal tissues using p< 0.05 and |log2FC| > 1.5 as cutoffs (Additional file 2, Table S2). These DEGs belong to 6 subclasses (LIM, ANTP, CUT, POU, PRD, TALE, Figure 2). Among 12 LIM homeobox genes (LHXs), 6 LHXs (LHX1, LHX2, LHX3, LHX5, LHX9, LMX1B) were indicated as upregulated, while ISL1 and LMX1A were found to be downregulated in HNSCC patients. For ISL2, LHX4, LHX6, and LHX8, no significant changes were observed (Figure 3A). The correlation of their expression between each other was calculated using Pearson’s correlation analysis and shown in Figure 3B. The gene with the most co-expression genes was LHX1, which was positively correlated with LHX5, LMX1B, LHX3, and negatively correlated with LHX2, ISL1, and LMX1A (Pearson’s correlation coefficient = 0.24, 0.23, 0.11, -0.26, -0.18, -1, respectively).
Next, we validated the altered expression of 8 LHXs in tumor and corresponding normal tissue using paired t-test. Similarly, significant up-regulation of LHX1, LHX2, LHX3, LHX5, LHX9, LMX1B and downregulated of ISL1 and LMX1A were examined in HNSCC tissues (all p<0.0001, Figure 3C). Furthermore, the ability to distinguish between normal and tumor samples was tested using the ROC curve. As shown in Figure 3D, all the selected LHXs showed favorable accuracy (LHX5: AUC = 0.956, LHX1: AUC = 0.919, LHX2: AUC = 0.913, LHX9: AUC = 0.815, LHX3: AUC = 0.757, LMX1B: AUC = 0.776, LMX1A: AUC = 0.750, ISL1: AUC = 0.735).
Association Between LHXs Expression and Clinicopathological Features
To understand the clinical significance of the selected LHXs, we analyzed the relationships between clinical variables and LHXs levels in patients with HNSCC. As shown in Figure 4, the expression of LHX1 and LHX2 was correlated with clinical stages (p= 0.041, p= 0.019, respectively) and patients with more advanced clinical stages have higher LHX2 expression levels (Figure 4A). As for LHX2 and LMX1A, high mRNA expression was correlated with advanced tumor grade (p= 7.23e-5, p= 0.01, respectively), in contrast to LHX1 which was highly expressed in G1 tumors (p= 0.015) (Figure 4B). Figure 4C displayed the significant correlations between LHX1, LHX2, LHX5 expression and primary tumor sites (all p< 0.0001). The highest expression of LHX1, LHX2, and LHX5 was evaluated in the oral cavity, pharynx and larynx, respectively. For HPV status, Figure 4D showed a remarkable higher expression for LHX1 and LHX5 in HPV-negative patients (p = 8.98e-13 and 1.57e-5). Conversely, LHX2 was found to express highly in HPV-positive patients(p=1.22e-10).
Prognostic Values of LHXs Expression in Patients With HNSCC
Figure 5 showed the association between patient survival (OS, DSS, DFS) and LHXs expression levels using the KM curve, LHX1 tended to have significantly adverse prognosis of OS, DSS, and DFS when expressed highly (OS and DSS, p < 0.0001; DFS, p = 0.0055). High expression of ISL1 was correlated with good prognosis for OS and DSS (p = 0.017 and 0.011, respectively), but not for DFS (p = 0.4). High LHX2 expression was associated with favorable OS and DFS (p = 0.048 and 0.02, respectively), but not DSS (p = 0.11). LHX5, as well as LMX1B, was correlated with unfavorable OS and DSS (LHX5: OS, p = 0.0053; DSS: p = 0.0014; LMX1B: OS, p = 0.00016; DSS, p = 0.0087), and showed no statistical association with DFS (LHX5: p = 0.31; LMX1B: p = 0.16). High level of LHX9 was related to good DSS and DFS (p = 0.032 and 0.038, respectively), but not OS (p = 0.17). On the contrary, high expression level of LMX1A was statistically associated with adverse prognosis in patients with HNSCC (DSS, p = 0.089; DFS, p = 0.03), but independent of OS (p = 0.21). There is no significant relationship between the expression of LHX3 and OS, DSS or DFS in patients with HNSCC.
The findings were further verified using Cox regression analysis. As shown in Figure 6, LHX1, LHX5, and LMX1B were risk factors for OS and DSS (OS: LHX1: HR, 2.15; 95% CI, 1.55-2.98; LHX5: HR, 1.58; 95% CI, 1.16-2.15; LMX1B: HR ,1.92; 95% CI, 1.36-2.71; DSS: LHX1: HR, 2.56; 95% CI, 1.73-3.78; LHX5: HR, 1.75; 95% CI, 1.24-2.57; LMX1B: HR, 1.75; 95% CI, 1.15-2.66; all p < 0.05), while ISL1 and LHX2 were related to favorable OS (ISL1: HR, 0.51; 95% CI, 0.3-0.86; LHX2: HR, 0.68; 95% CI, 0.46-1; all p < 0.05). ISL1 and LHX9 were predictors of favorable DSS (HR, 0.42; 95% CI, 0.21-0.84; HR, 0.68; 95% CI, 0.48-0.97; respectively, all p < 0.05). In terms of DFS, LHX1 and LMX1A were indicated as unfavorable predictors (HR, 4.7; 95% CI, 1.41-15.67; HR, 2.31; 95% CI, 1.06-5.05; respectively, all p < 0.05), while LHX2 and LHX9 were indicated as favorable factors (HR, 0.13; 95% CI, 0.02-0.97; HR, 0.41; 95% CI, 0.17-0.98; respectively, all p < 0.05).
LHXs Correlated Genes Are Involved in Numerous Biological Processes
After indicated selected LHXs as prognosis factors for patient survival, we conducted the GSEA analysis to assess these genes at a functional level (Additional file 2, Table 3). We found that high-LHX1 expression exhibited enrichment of genes associated with environmental information processing pathway (ECM-receptor interaction and ) and immune system pathway (RIG-I-like receptor signaling pathway, cytosolic DNA-sensing pathway ). Gene sets that are involved in metabolism processes like metabolism of xenobiotics by cytochrome P450 were associated with low expression of LHX1 and LMX1B. As for LHX5, GSEA analysis revealed that pathways like focal adhesion and ribosome were up-regulated and pathways related to immune response (T cell receptor signaling pathway, primary immunodeficiency, and Natural killer cell mediated cytotoxicity) were down-regulated. Pathways implicated in cancer invasion and metastases like focal adhension and tight junction were activated in the high-expression groups of LHX5, LMX1A, LMX1B and low-expression groups of LHX2 and LHX9. The enrichment of post-transcriptional control related pathways like ribosome, spliceosome, proteasome, and RNA degradation was correlated to high LHX3 and LHX5 expression while low ISL1, LHX2 and LMX1A expression. Moreover, LHXs have a close relationship with immune-related pathways, which were enriched among the low-expression groups of ISL1 (cytosolic DNA sensing pathway, antigen processing and presentation, RIG-I-like receptor signaling pathway), LHX2 (cytosolic DNA sensing pathway), LHX9(intestinal immune network for IgA production, T-cell receptor signaling pathway), and LMX1B (primary immunodeficiency), Figure 7A.
Next, differential gene sets associated with high and low LHXs expression phenotypes with regard to hallmark gene sets (h.all.v7.4) were analyzed (Additional file 2, Table S4). High LHX1 expression was related to the process of IL6 JAK STAT3 signaling, TGF beta signaling and angiogenesis. High expression of LHX3 was a connection with G2M checkpoint, MYC family targets, DNA repair, TGF alpha response, and E2F targets. LHX5 and LMX1B high expression was associated with angiogenesis and epithelial-mesenchymal transition (EMT). Low ISL1 expression corresponded with the process of reactive oxygen species (ROS) pathway, TGF beta signaling and MYC targets v2. Low LHX2 expression corresponded with angiogenesis, epithelial-mesenchymal transition, P53 pathway, TNF signaling via NF-kB, stress reaction pathways like unfolded protein response and UV response, immune response pathways like inflammatory response, interferon-alpha response, and interferon-gamma response. Low LHX9 expression corresponded with the process of IL2 STAT5 signaling, KRAS signaling and immune response pathways like allograft rejection and inflammatory response, Figure 7B.
Immune Infiltration Is associated with LHXs Expression
Figure 8 displayed the CIBERSORT analysis of the relationship between tumor-infiltration immune cells and LHXs expression. ISL1 expression level had negative associations with the ratio of resting memory CD4+ T cells, M0 macrophages, M1 macrophages, and positive relationships with the ratio of plasma cells. The expression levels of LHX1 and LMX1B showed positive correlations with the proportion of resting memory CD4+ T cells, resting NK cells, CD4+ T cells, M0 macrophages, as well as M1 macrophages, and negative association with the ratio of plasma cells and CD8+ T cells. LHX2 expression was positively correlated with the infiltration of plasma and CD4+ T cells, but negatively correlated with M0 macrophages and activated dendritic cells. A higher proportion of resting NK cells and M1 macrophages, while a lower proportion of activated dendritic cells was observed in the LHX3 high-expression group. LHX5, as well as LMX1A, had high M0 macrophages frequency and low plasma cells, CD4+ T cells, M1 macrophages frequencies in their high-expression group. As for LHX9, the CIBERSORT output revealed that the M1 macrophages and activated dendritic cells proportion were significantly higher in LHX9 high-expression group.
The tumor-infiltration immune cells according to LHXs expression evaluated in Timer database was shown in Figure 9. The increased immune cells infiltration in total HNSCC patients was found to correlated with the expression of ISL1 (B cell: cor = 0.209, p = 4.07e-6; CD8+ T cell: cor = 0.097, p = 3.51e-2; CD4+ T cell: cor = 0.143, p = 1.71e-3), LHX2 (B cell: cor = 0.181, p = 6.96e-5; CD8+T cell: cor = 0.157, p = 6.02e-4; CD4+ T cell: cor = 0.202, p = 8.18e-6; Macrophage: cor = 0.111, p = 1.43e-2; Neutrophil: cor = 0.138, p = 2.47e-3; Dendritic cell: cor = 0.191, p = 2.41e-5), and LMX1A (B cell: cor = 0.158, p = 5.46e-4; CD8+ T cell: cor = 0.137, p = 2.78e-3; CD4+ T cell: cor = 0.253, p = 1.99e-8; Macrophage: cor = 0.31, p = 3.30e-12; Neutrophil: cor = 0.11, p = 1.59e-2; Dendritic cell: cor = 0.198, p = 1.16e-5). The decreased immune cells infiltration was evaluated to be related to the expression of LHX5 (B cell: cor = -0.139, p = 2.36e-3; CD8+ T cell: cor = -0.19, p = 3.15e-5; Neutrophil: cor = -0.119, p = 9.35e-3; dendritic cell: cor = -0.108, p = 1.72e-2) and LHX9 (B cells: cor = -0.119, p = 1.77e-2; Macrophages: cor = -0.188, p = 3.32e-5; Dendritic cell: cor = -0.079, p = 8.41e-2). LHX1 was positively related to CD4+ T cell, neutrophil, and dendritic cell infiltration (cor = 0.117, p = 1.04e-2; cor = 0.177, p = 9.42e-5; cor = 0.094, p = 3.81e-2; respectively), but negatively related to B cell and CD8+ T cell infiltration (cor = -0.194, p = 2.04e-5; cor = -0.159, p = 5.14e-4; respectively). LMX1B was examined with positive association a decrease of B cell (cor = -0.145, p = 1.52e-3) infiltration, as well as an increase of Macrophage (cor = 0.098, p = 3.08e-2) infiltration.
Relationship between DNA Methylation and LHXs Dysregulation
To understand the epigenetic regulatory mechanisms underlying the abnormal expression of LHXs in HNSCC, we explored the methylation level of LHXs promoters in the UALCAN database. The significant hyper-methylation of ISL1 (p = 1.624e-12), LHX2 (p = 5.384e-8), LHX5 (p < 1e-12), LHX9 (p < 1e-12), LMX1A (p < 1e-12), and LMX1B (p < 5.754e-4) promoters and hypo- methylation of LHX3 (p = 5.384e-8) promoter was displayed in Figure 10A. There is no statistically significant change in the methylation level of LHX1 promoters between HNSCC tissues and normal samples (p = 0.145).
We further explored the relationship between DNA methylation and the mRNA expression levels of LHXs that are down-regulated in HNSCC tissues using MEXPRESS. As shown in Figure 10B, low ISL1 expression was found to be significantly associated with the hyper-methylation for all promoter regions probes in HNSCC tissues (Pearson correlation coefficients: −0.329 ~ −0.537, all p < 0.001). As for LMX1A, we found the 3 significant methylation sites (cg09248345, cg04351049, and cg09069138) in the DNA sequences of LMX1A that were negatively associated with its expression levels (Figure S1).