Evaluation of the HOXA11 level in patients with lung squamous cancer and insights into potential molecular pathways via bioinformatics analysis

Background This study was carried out to discover the underlying role that HOXA11 plays in lung squamous cancer (LUSC) and uncover the potential corresponding molecular mechanisms and functions of HOXA11-related genes. Methods Twenty-three clinical paired LUSC and non-LUSC samples were utilized to examine the level of HOXA11 using quantitative real-time polymerase chain reaction (qRT-PCR). The clinical significance of HOXA11 was systematically analyzed based on 475 LUSC and 18 non-cancerous adjacent tissues from The Cancer Genome Atlas (TCGA) database. A total of 102 LUSC tissues and 121 non-cancerous tissues were available from Oncomine to explore the expressing profiles of HOXA11 in LUSC. A meta-analysis was carried out to further assess the differential expression of HOXA11 in LUSC, including in-house qRT-PCR data, expressing data extracted from TCGA and Oncomine databases. Moreover, the enrichment analysis and potential pathway annotations of HOXA11 in LUSC were accomplished via Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The expression of hub genes and according correlations with HOXA11 were assessed to further explore the biological role of HOXA11 in LUSC. Results HOXA11 expression in LUSC had a tendency to be upregulated in comparison to adjacent non-cancerous tissues by qRT-PCR. TCGA data displayed that HOXA11 was remarkably over-expressed in LUSC compared with that in non-LUSC samples, and the area under curves (AUC) was 0.955 (P < 0.001). A total of 1523 co-expressed genes were sifted for further analysis. The most significant term enriched in the KEGG pathway was focal adhesion. Among the six hub genes of HOXA11, including PARVA, ILK, COL4A1, COL4A2, ITGB1, and ITGA5, five (with the exception of COL4A1) were significantly decreased compared with the normal lung tissues. Moreover, the expression of ILK was negatively related to HOXA11 (r = − 0.141, P = 0.002). Conclusion High HOXA11 expression may lead to carcinogenesis and the development of LUSC. Furthermore, co-expressed genes might affect the prognosis of LUSC.

Background Non-small cell lung cancer (NSCLC) has caused the most frequently cancer-related deaths among all types of malignancy in humans worldwide, accompanied by a high incidence [1][2][3]. NSCLC is responsible for the majority of the primary lung cancer cases, including large cell carcinoma, lung adenocarcinoma (LUAD), and lung squamous cancer (LUSC). Among the three histological subtypes, LUSC is the most common type in developing countries [4][5][6][7], and patients with lung cancer are still facing a low overall 5-year survival rate [8,9]. Molecular targeted therapy has achieved curative efficacy in clinical in LUSC. For example, EGFR targeted therapy had a modest effect in advanced LUSC patients [10]. However, the number of applicable patients is limited [11]. Thus, there is an urgent need to identify more underlying highperformance targets in LUSC.
Cumulative evidence has demonstrated that HOX genes, which belong to the large family of homeodomain genes, work to regulate growth processes, such as organogenesis and body patterning [12,13]. Humans have HOX genes in four clusters (HOXA, HOXB, HOXC, and HOXD) [14,15]. Cumulative studies have reported that HOX genes are expressed in healthy human lungs and play a crucial role in their development [12]. The HOXA11 cluster is located on chromosome 7p15-7p14.2, and 12 genes are involved in the cluster, including EVX1 and 11 HOX genes [16,17]. A number of studies have been carried out to define the function of HOXA genes in malignant cancers. An increasing number of reports suggest that HOXA11 has been implicated in several malignant tumors, such as gastric cancer [18], renal cell carcinoma [19], NSCLC and lung adenocarcinoma [17,20], and breast cancer [21]. Thus far, the expression level of HOXA11 and its potential mechanisms in LUSC have not been clarified.
In the present study, we attempted to identify the association between clinical parameters and HOXA11 expression to gain a comprehensive understanding of the role of HOXA11 in LUSCs. The mechanisms of HOXA11 co-expressed genes were mined by bioinformatics analysis.

Selection of clinical LUSC tissue samples
Clinical samples were collected from 23 LUSC patients who had been pathologically identified at the Department of Pathology, First Affiliated Hospital of the Guangxi Medical University (Nanning, Guangxi, China) , from January 2012 to February 2014. The clinicopathological features of the patients are shown in Table 1.
The Ethical Committee of the First Affiliated Hospital of Guangxi Medical University approved the present research. All participating clinical doctors and patients signed written informed consents.

Total RNA isolation
Per the manufacturer's instructions, we extracted total RNA with the miRNeasy FFPE Tissue Kit (QIAGEN, Shanghai, China). In addition, we detected the purity and concentration of total RNA using NanoDrop 2000 (ThermoScientific, USA).

Data mining and analyzing
All clinicopathological parameters related to LUSC and mRNA (level 3) expression in LUSC were carefully downloaded from the TCGA Data Portal website (http:// cancergenome.nih.gov). Based on the HOXA11 expression in LUSC, GraphPad Prism was applied to obtain the scatter diagram. In addition, SPSS was carried out to acquire receiver operating characteristic curves (ROCs) as well as overall survival (OS) and disease-free survival (DFS) curves. Meanwhile, the available data of HOXA11 expression in LUSC was mined in Oncomine (https:// www.oncomine.org). Furthermore, we collected in-house qRT-PCR and data from public databases to gain insight into the differential expression of HOXA11 in LUSC using an integrative meta-analysis. The standard mean difference (SMD) was pooled from all studies to determine the expression level of HOXA11 in LUSC.

Enrichment analysis and pathway annotation
Gathered genes were analyzed using bioinformatics. The enrichment of functions and signaling pathways of the target genes were analyzed using The Database for Annotation, Visualization and Integrated Discovery v6.8 (DAVID), FunRich, and Cytoscape. The String database (http://www.string-db.org) was applied to construct the protein-protein interaction (PPI) network for the hub gene identification. Moreover, hub genes were selected to obtain their expression and correlation with HOXA11 in LUSC. The immunohistochemistry results of the six hub genes in LUSC were retrieved from the Human Protein Atlas (HPA) database.

Statistical analysis
Statistical analysis was conducted using SPSS 22.0. All data are shown as mean ± standard deviation (SD). An independent samples' t test was adopted to examine the differences between cancer tissues and normal lung tissues. A one-way analysis of variance (ANOVA) was employed for analyzing differences of HOXA11 expression in various pathological gradings from inhouse qRT-PCR data, as well as terms of M category, N category, race, and statuses of recurrence from the TCGA database. The relationships between the coexpressed genes and HOXA11 were assessed using the Pearson rank correlation, and the AUC was counted.
To achieve an in-depth understanding of the prognostic value of HOXA11, we also used Kaplan-Meier curves to determine the survival time, including the OS and DFS. A two-sided P value < 0.05 was considered statistically significant.

HOXA11 expression and clinicopathological features in LUSC using qRT-PCR
The expression and clinicopathological features of HOXA11 in LUSC are displayed in Table 1. There was no significant correlation between HOXA11 expression  and all clinical parameters. However, the expression of HOXA11 was upregulated in LUSC compared to in non-cancerous tissues, and the AUC of the TNM stage was 0.831 (P = 0.008) (Fig. 1a, b).

Bioinformatics analysis
All analysis was based on the number of 1523 coexpression genes. The Gene Ontology enrichment analysis comprised three categories: a biological process (BP), a molecular function (MF), and a cellular component (CC). The most valuable 10 pathways of each category are presented in Fig. 6a, c, and 6e, including different kinds of functional relationship graphs (Fig. 6b, d, and f ). For the KEGG pathways, the 10 most significant pathways are shown in Table 4. The PPI network is displayed (Fig. 7); three pairs of hub genes with the highest combined scores (PARVA, ILK, COL4A1, COL4A2, ITGB1, and ITGA5) were collected from the PPI network (Table 5).

Discussion
HOX genes may play a central role in regulating gene expression, differentiation, and receptor signaling. Members of HOX family can also encode DNA-binding transcription factors [30,31]. A growing body of studies observed that HOXA11 expression was downregulated in different tumor types. In glioblastoma, the decreased level of HOXA11 was confirmed as a significant marker of poorer prognosis. [32]. Moreover, Bai et al. suggested that the down-expressed HOXA11 gene may play an essential role in carcinogenesis by promoting gastric cancer development, which may be helpful to forecast the malignant behaviors of gastric cancer [33]. HOXA11 has been notably over-expressed in epithelial ovarian cancers [34]. Data from the current study, TCGA, Hou's study [25], along with Garber's study [26], all illustrated that the corresponding mRNA expression of HOXA11 is significantly over-represented in LUSC compared with adjacent tissues. The enhanced expression of HOXA11 might be a diagnostic target to use for distinguishing LUSC from healthy controls based on the ROC (AUC = 0.955, P < 0.001) from TCGA. Unfortunately, the expression detected by qRT-PCR revealed no significance. However, HOXA11 tends to be upregulated in LUSC compared to normal lung tissues. Interestingly, Talbot reported that the levels of HOXA11 were reduced in patients with LUSC. In light of these previous studies, HOXA11 expression seems to be higher in LUSC than in healthy controls. More LUSC specimens need to be collected to study the expression levels of HOXA11 in LUSC in the future. Among the co-expressed genes of HOXA11, a total of three genes, HOXA10, HOXA13, and HOXC10, were finally intersected, which all belong to the HOX gene family. HOXA10 is involved in gene expression, regulation, morphogenesis, and differentiation in ovarian carcinoma. HOXA10 and HOXA11 might be associated with primary tumors and specific histological subtypes [34]. HOXA13 regulates gene expression and differentiation. HOXA10, HOXA11, and HOXA13 might be useful targets to further mine the molecular pathogenesis of HOXA11 in early-stage lung adenocarcinoma [35]. Moreover, HOXC10 may act as the accelerator for original activation. In previous studies, HOXC10 has been reported to be associated with the increased invasion of malignancies [36,37]. These three co-expressed genes and HOXA11 might play several pivotal roles in LUSC, such as the subtype differentiation of lung cancer, the regulation of LUSC progression, and the development of efficient therapeutic strategies.
The pathway of focal adhesion has been extensively reported to be the most significant pathway in various diseases and signaling pathways as well. It regulates diverse cellular functions that were once activated by focal adhesion kinase (FAK), including adhesion, proliferation, migration, and survival [38][39][40]. A pattern Fig. 6 Top 10 significant pathways and GO enrichment analysis. a Graph of the 10 most significant pathways of BP category. b Enrichment analysis of BP, each node means one different function and more significant ones are filled in with a deeper color. c Top 10 significant terms in the CC category. d Enrichment analysis of the CC category; each node means one different function, and more significant ones are filled in with a deeper color. e Ten most valuable annotations of the MF category. f Enrichment analysis of the CC category; each node means one different function, and more significant ones are filled in with a deeper color Table 4 Top 10 significant KEGG pathway annotations   Term   Genes   Count  P value   Focal adhesion  CAV2, CAV1, TNC, COL3A1, ITGA11, PTEN, ITGB1, DOCK1, LAMB2, COL27A1, ITGAV, ILK, COL6A2, COL6A1, PDGFC, PDGFD, LAMB1, THBS1, THBS2,  AKT3, FN1, PRKCA, COL4A2, COL4A1, ROCK2, MET, ITGA1, ITGA4, BAD, FLNC, COL5A2, COL5A1, FLNA, LAMA2, VEGFC, LAMA4, ITGA5, COL1A2,  PDGFRA, PDGFRB, COL1A1, LAMC1, MYLK, PARVA   44   1.21E−12   ECM-receptor interaction  TNC, COL3A1, ITGA11, ITGB1, LAMB2, CD44, ITGAV, COL27A1, COL6A2, COL6A1, THBS1, LAMB1, THBS2, FN1, COL4A2, COL4A1, HSPG2, ITGA1, ITGA4,  COL5A2, COL5A1, LAMA2, LAMA4, ITGA5,   of enhanced expression of FAK in lung carcinomas has been reported, which is related to nodal involvement and the deterioration of advanced disease stages [41,42]. Therefore, FAK protein expression may help in predicting the aggressive behavior of LUSC. Meanwhile, FAK might be pursued as a promising therapeutic target for LUSC. Six hub genes were analyzed, PARVA, ILK, COL4A1, COL4A2, ITGB1, and TIGA5. Except COL4A1 gene, the rest of the hub genes were all significantly decreased in LUSC compared to adjacent healthy controls. COL4A1 and COL4A2 both belong to the type IV collagen gene family. Both can act as inhibitors of angiogenesis and tumor growth. In LUSC, they can control vascular invasion and inhibit the size of tumors. There is more reason to believe that the upregulated HOXA11 gene plays a carcinogenic role in LUSC. HOXA11 might depress COL4A1 and COL4A2 expression levels in LUSC. ITGB1 has been recognized in the processing of metastatic diffusion of tumor cells, and ITGA5 may promote tumor invasion. In lung cancer, higher expression of ITGA5 may be correlated with a shorter survival time. Meanwhile, previous studies have indicated that modulating the ILK signaling pathway by PARVA made it more vulnerable to metastasis in lung adenocarcinoma [43]. These hub genes might also similarly interact with each other via various signaling pathways in LUSC. Thus, we speculate that LUSC patients' prognosis might be better and their survival time might be longer than for patients with other subtypes of lung cancer. The mechanisms and functions between hub genes and HOXA11 in LUSC remain elusive and need to be validated in the future.

Conclusion
The findings of the present study indicate that HOXA11 high expression might lead to the occurrence and development of LUSC. Meanwhile, co-expressed HOXA11 genes may influence the prognosis of LUSC, and the gene ILK may have the complete reverse functions in LUSC compared with HOXA11. Hub genes need to be further analyzed to ensure their mechanisms and functions in LUSC.  8 Hub genes' expression in LUSC and correlations with HOXA11. a ILK was lower in LUSC tissues than in non-cancerous tissues (P < 0.001). b The gene PARVA was significantly overexpressed in normal tissues (P < 0.001). c The levels of COL4A1 in different tissues showed no significance (P = 0.061). d The hub gene ITGB1 revealed higher levels in normal tissues (P < 0.001). e ITGA5 was significantly decreased in LUSC tissues (P < 0.001). f COL4A2 upregulated in non-LUSC tissues (P < 0.001). g ILK and HOXA11 showed a negative correlation (r = − 0.141, P = 0.002).