Abnormal expression of HOXD11 promotes the malignant behavior of glioma cells and leads to poor prognosis of glioma patients

Data collection

GEPIA is an integrated database platform, which is characterized by the fusion of gene sequencing data of tumor samples from TCGA database and gene sequencing data of normal control group sample tissues from The Genotype-Tissue Expression (GTEx) (Tang et al., 2017). This advantage of GEPIA makes up for the shortage of normal control samples in TCGA database. The expression levels of HOXD11 in various tumors and tissues, including 163 glioblastoma, 518 low grade glioma and 207 normal brain tissue samples, were obtained from the GEPIA database. To further determine the expression level changes of HOXD11 in glioma, we obtained gene expression data in the form of FPKM (Fragments Per Kilobase per Million) from a glioma tissue microarray (GSE4290) and a glioma cell microarray (GSE15824) from GEO database https://www.ncbi.nlm.nih.gov/geo/ for further analysis. GSE4290 contains 23 non-tumor brain tissue samples and 77 glioma tissue samples (Sun et al., 2006). GSE15824 contains three human-derived astrocyte (HA) cell line samples and two LN319, two LN229, two BS149 and two LN018 glioblastoma cell line samples (Grzmil et al., 2011).

CGGA (http://www.cgga.org.cn/) is a public database focused on glioma research, which contains many types of data and corresponding clinical information (Guan et al., 2018). This part of the data contains a variety of clinical information of the patients and was used to further depth analyze the relationship between changes in HOXD11 expression level and glioma prognosis. We examined two RNA-seq datasets containing 693 and 325 samples. After excluding the data of patients with incomplete clinical data, we obtained 749 glioma samples for further analysis of various types.

In order to verify the effect of HOXD11 on the prognosis and diagnostic value of glioma, we searched and obtained one RNA-seq dataset containing 666 glioma tissues in the TCGA database (https://portal.gdc.cancer.gov/) and obtained clinical information of the corresponding patients.

Nine glioma tissue samples and five non-tumor brain tissue samples were collected and stored in a liquid nitrogen environment during surgery and then transferred to a -−80 °C freezer until use. Clinical samples were obtained after receiving written informed consent from patients in the clinical department of Henan Provincial People’s Hospital. This study has been approved by the Institutional Review Board of Zhengzhou University.

GSEA analysis of HOXD11

GSEA is widely used as a bioinformatics analysis tool for gene function annotation and analysis. The RNA-seq data obtained from the CGGA database was batch-corrected and normalized by limma packages and then classified into “group H” high expression group or “group L” low expression group based on the median expression level of HOXD11. Functional enrichment analysis was performed using the GSEA 4.0.2 jar software. The number of permutations was set to 1000 times, and the gene set database was set to the Kyoto Encyclopedia of Genes and Genomes (KEGG) cell signaling pathway.

Cell culture

Human glioma cell lines (U251, LN229 and T98) and human-derived astrocyte (HA) were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). All cells were grown in incubators at 37 °C and 5% carbon dioxide and cultured using DMEM medium (Thermo Fisher Scientific, USA) plus 10% FBS (Thermo Fisher Scientific, USA).

RNA extraction and quantitative reverse transcription polymerase chain reaction (RT-qPCR) analysis

9 glioma tissue and 5 peritumoral tissue brain tissue samples were used to detect the expression level of HOXD11 in tissues by RT-qPCR. Three glioma cell lines (T98, U251, LN229) and one human astrocyte were used to detect the expression level of HOXD11 in glioma cells. First, total RNA was extracted from glioma and normal brain tissue samples using Tri-Reagent (Sigma, USA). The quality and quantity assessments of total RNA were determined using NanoDrop One spectrophotometer (Thermo Fisher Scientific, USA), evaluating 260/280 nm absorbance values. Thereafter, the cDNA was reverse transcribed from the total RNA using the Transcriptor First Strand cDNA Synthesis kit (Roche, USA). RT-qPCR was performed following the guidelines for FastStart Universal SYBR Green Master (ROX) (Roche, Germany). Results were quantified using QuantStudio software (Thermo Fisher Scientific, USA) following the manufacturer’s instructions. GADPH was used as an internal reference. The primer sequences used for HOXD11 were 5′-GAACGACTTTGACGAGTGCG-3′(F) and 5′-ACGGTTGGGAAAGGAACGAA-3′(R). The expression level of HOXD11 was measured by the “2–ΔΔCT” method. Statistical differences between the two groups were analyzed by unpaired t-test, and results were considered statistically significant when the p value was < 0.05.

Cell transfection

The small interfere RNA (siRNA) that specifically targeting HOXD11 was purchased from Genepharma (Shanghai, China). The siRNA sequence was sense (5′-GACUUCAACUCUCUCGGAUTT-3′) and antisense(5′-AUCCGAGAGUUGAAGUCTT-3′). The negative control (NC) RNAi sequence was sense(5′-UUCUCCGAACGUGUCACGUTT-3′) and antisense, (5′-ACGUGACACGUUCGGAGAATT-3′). Cells were transfected using Lipofectamine 2000 reagent (Invitrogen, Carlsbad, CA, USA). Detection of transfection efficiency by RT-qPCR at 48 h after transfection.

Cell proliferation assay

To detect the effect of HOXD11 on the proliferation ability of U251 cells, MTT experiments were conducted. The transfected cells were transferred into 96-well plates (1 × 10³ cells/well). Cell proliferation ability was measured at 12 h, 24 h, 48 h and 96 h after transfection. The specific method was to add 20 ul of MTT (5mg/ml) to each well for culture. The liquid was discarded after incubation for 4 h and 150 ul DMSO was added to each well. After 15 min of shaking, the optical density (OD) of 490 nm was detected using microplate reader.

Colony formation assay

In colony formation assay, U251 cells were planted in six-well plates at a density of about 300 cells/well. After 10 days of culture, the cells were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Colonies larger than one mm in diameter were counted using ImageJ Software (v1.8.0) after photographing.

Wound-healing assay

The exponential growth cells were seeded on the 6-well plate with a fusion degree of 90%. A straight linear scratch was made on the center of the plate with a pipette (200 uL). PBS was washed gently, and serum-free medium was replaced for culture. After that, the scratch width was observed and photographed under the microscope for 24 h and 48 h respectively, then was measured by ImageJ Software (v1.8.0).

Immunofluorescence staining

Firstly, treated cells were fixed with 4% paraformaldehyde for 30 min, and then 0.3% Triton X-100 was used to penetrate for 20 min at room temperature. After blocking with 5% BSA, primary antibody Ki67 (1:1000, Abcam, UK) was added and placed in a 4 °C refrigerator overnight. After that, washed 3 times with PBS for 5 min each time. Secondary antibody was added and incubated for one hour at room temperature in the dark. After washing with PBS, DAPI was added for nucleation. The photographs were finally observed using a fluorescence microscope. After that, the total number of cells and the number of KI67 positive cells were counted by ImageJ Software (v1.8.0).

Flow cytometry

U251 cells were digested and centrifuged to collect the cells after transfection. Precooled 75% ethanol was added, fixed overnight at −20 °C Each tube was added with 500 uL propidium iodide (PI) dye solution, stained for 30 min in dark, and filtered through a filter. Cell cycle distribution was analyzed by flow cytometry (Canto plus, BD, USA).

Western blotting

Total protein of transfected U251 cells was extracted using RIPA buffer with protease inhibitor. Separated by 10% SDS-PAGE gel electrophoresis, transferred to polypropylene fluoride membrane (PVDF) treated with methanol for 2 h, transferred to a shaker, sealed with TBS-T blocking solution containing 5% skimmed milk powder for 1 h. Then, the primary antibodies GAPDH (1:5000, Abcam, UK), CCNB1 (1:500, Abcam, UK), CDK1 (1:500, Abcam, UK) and Cdc25C (1:500, Abcam, UK) were added to incubate overnight at 4 °C. Washed by TBS-T 5 times for 5 min each time, secondary antibody (1:5000) was added to incubate for 1 h at room temperature, Washed by TBS-T 5 times for 5 min each time, and finally ECL color developing solution was added for color development. The developed bands were subjected to gray level analysis by ImageJ Software (v1.8.0), The difference in protein expression levels between the HOXD11-NC group and the HOXD11-siRNA group is represented by the gray value of the corresponding protein bands.

Statistical analysis

R software (v.3.6.1 version) and GraphPad Prism (v8.0.2) was used to process data and perform statistical analysis. Unpaired t-test was used for comparison between the two groups. The expression levels of HOXD11 were examined in glioma and non-tumor brain tissue samples by the Wilcox test method. The relationship between the Overall Survival rate of patients and the expression level of HOXD11 was analyzed by the Kaplan–Meier method and Cox regression analysis, and the survival curve was drawn accordingly. The relationship between clinically relevant patient information and the expression level of HOXD11 was analyzed using the Wilcox and Kruskal tests.

Abnormally high expression of HOXD11 in various tumors including gliomas

Through GEPIA database, we found that HOXD11 was abnormally highly expressed in a variety of tumors, such as ESCA, GBM, HNSC and HUSC, especially in GBM, as shown in Fig. 1A. To confirm the expression level of HOXD11 in gliomas relative to normal controls, data from more sources were analyzed. By analyzing the data of glioma tissues and cell lines microarrays obtained from GEO database, it was found that HOXD11 was highly expressed in both glioma tissues and cell lines, as shown in Figs. 1B and 1C. To validate the above analysis results, we performed RT-qPCR at tissue and cell levels, respectively. The results showed that HOXD11 was highly expressed in glioma tissues and cell lines compared with normal brain tissues and human astrocytes, as shown in Figs. 1D and 1E. Among them, the relevant clinical information of 9 glioma tissue samples was shown in Table S1.

Correlation between the expression level of HOXD11 and the clinical and prognosis of patients with glioma

To clarify the relationship between high expression of HOXD11 and the clinical and prognosis of glioma patients, we conducted further studies using glioma sequencing data and relevant clinical information from CGGA database. A total of 749 glioma samples with complete clinical information were obtained from the CGGA database, which contains a variety of clinical information such as Primary-Recurrent-Secondary (PRS) type, age (age at diagnosis), gender, histological type, and chemotherapy status, as well as 1p19q codeletion status and isocitrate IDH mutation status. More detailed information of clinical characteristics is shown in Table S2.

The Wilcox and Kruskal tests were used to analyze the relationship between clinically relevant patient information and the expression level of HOXD11. The expression level of HOXD11 was significantly correlated with PRS type, WHO grade, chemotherapy status, age, 1p19q codeletion data, IDH mutation status and histological type (Figs. 2A–2G). The expression level of HOXD11 increases as the WHO grade and age of the glioma increases (p < 0.001). As for the IDH mutation status, there was a higher level of gene expression in wildtype than in the mutant (p < 0.001). Moreover, for 1p19q codeletion status, there were lower gene expression levels in patients with a codeletion of 1p19q than in patients with a non-codeletion of 1p19q (p < 0.001).

Survival outcomes of HOXD11 in patients with glioma

The relationship between overall survival and the expression level of HOXD11 in patients with gliomas was explored by Kaplan–Meier survival analysis. The 749 patients with gliomas included in the study were divided into two groups: “group H” and “group L,” based on the median expression level of HOXD11. As shown in Fig. 3A, the Overall Survival (OS) was significantly lower in patients with gliomas having high expression levels of HOXD11 (group H) compared to those with low levels (group L; p < 0.001) in CGGA-sequence, which indicated poor prognosis.

To improve the credibility of the data results, we further obtained 666 glioma tissues in the TCGA database. The Kaplan–Meier method was used to verify that the increase of HOXD11 expression level can indeed reduce the overall survival time of glioma patients (Fig. 3B). To clarify the impact of HOXD11 expression level on the prognosis of gliomas of different molecular subtypes, we divided the samples into three groups (IDH mutation with 1p19q codeletion, IDH mutation without 1p19q codeletion and IDH wild type) for survival analysis. The results showed that the prognosis of the group with high HOXD11 expression was significantly worse than that of the group with low HOXD11 expression among the three molecular subtypes, as shown in Figs. 3C–3E. These results suggest that high expression of HOXD11 may lead to poor prognosis of patients with glioma.

High expression of HOXD11 is an independent risk element in patients with gliomas

To determine whether HOXD11 could be used as an independent prognostic factor, we performed univariate and multivariate analyses for the date from CGGA-sequence. Univariate analysis was performed by applying Cox regression model (Fig. 4A). The results suggested that high expression levels of HOXD11 was significantly related to poor prognosis as a risk factor for glioma (p < 0.001; HR = 1.550; 95% CI [1.456–1.650]). The PRS type, histology, WHO grade, age and other factors showed similar results. Next, multivariate analysis was performed using the Cox regression model (Fig. 4B). High expression of HOXD11 (p < 0.001; HR = 1.156; 95% CI [1.072–1.246]), PRS type (p < 0.001; HR = 1.947; 95% CI [1.655–2.291]), and high WHO grade (p < 0.001; HR = 2.673; 95% CI [1.955–3.655]) were significantly associated with poor prognosis. Together, the univariate and multivariate analyses results suggest that high expression levels of HOXD11 may be an independent risk factor for poor prognosis.

High expression of HOXD11 has certain evaluation value for glioma prognosis

To clarify the diagnostic value of highly expressed HOXD11 for glioma prognosis, we plotted the receiver operating characteristic (ROC) curve based on the sequencing data of glioma patients from CGGA and TCGA databases. As shown in Figs. 4C–4D, the AUC values were greater than 0.700 in each dataset, indicating that high expression of HOXD11 had certain evaluation value for one-year, three-year and five-year survival rates. These results indicate that HOXD11 expression level could be related to prognosis in patients with gliomas; moreover, high expression of HOXD11 may be used for clinical prognosis evaluation of patients with glioma.

Knockout of HOXD11 could significantly inhibit the proliferation ability of U251 cells

By analyzing the sequencing data of glioma samples and matching clinical information from multiple sources, the results showed that HOXD11 could be an independent prognostic factor for glioma and had certain diagnostic value. To verify the effect of HOXD11 on malignant biological behavior of glioma cell lines, we conducted a series of phenotypic experiments. First, we used siRNA to interfere with the expression level of HOXD11 in glioma cell U251. Meanwhile, the transfection efficiency was detected by RT-qPCR (Fig. 5A). To determine the effect of HOXD11 on the proliferation ability of U251 cells, we performed MTT assay. The results showed that the proliferation ability of HOXD11 knockdown group(HOXD11-siRNA group) was significantly reduced compared with the control group(HOXD11-NC) (Fig. 5E). Meanwhile, HOXD11 knockdown also reduced the ability of cell colony formation (Fig. 5B). Meanwhile, immunofluorescence staining results showed that the expression level of cell proliferation marker KI67 in the knockdown group was significantly lower than that in the control group (Figs. 5C, 5D). In addition, similar results were obtained in the wound-healing assay, and the wound healing ability was also significantly reduced in the HOXD11 knockdown group (Fig. 5F). Therefore, we believe that highly expressed HOXD11 is involved in the pathological process of glioma by affecting the proliferation and migration of glioma cells.

Determination of HOXD11-related cellular signaling pathway by GSEA

To clarify how HOXD11 plays a role in the pathological process of glioma, we performed GSEA analysis to predict the cell signaling pathways that HOXD11 may participate in. As shown in Table 1 and Fig. 6A, HOXD11 may be involved in a variety of cancer-related cell signaling pathways such as cell cycle, DNA replication, ECM receptor interaction, and focal adhesion. These results reveal that HOXD11 may be involved in these cancer-related cell signaling pathway as an oncogene.

High expression of HOXD11 participates in the pathological process of glioma by regulating cell cycle

Through GSEA analysis, we predicted that highly expressed HOXD11 may participate in the pathological process of glioma through a variety of cancer-related cell signaling pathways.

We further chose cell cycle signaling pathway as an object to verify the prediction results. To explore the effect of HOXD11 on the cell cycle of glioma cell U251, we performed flow cytometry. The results showed that in the HOXD11 knockdown group, U251 cells were significantly blocked in G2 phase compared with the control group (Fig. 6B). To validate this result, we performed western blotting. The results showed that the expression levels of cell cycle-related proteins (CCNB1, CDK2 and Cdc25c) were significantly reduced in U251 cells with HOXD11 knockdown, confirming that HOXD11 is involved in regulating cell cycle progression (Fig. 6C).