DLX1 acts as a novel prognostic biomarker involved in immune cell infiltration and tumor progression in lung adenocarcinoma

TIMER database analysis

Tumor Immune Estimation Resource (TIMER, http://timer.cistrome.org/) is an online database for the comprehensive analysis of tumor-infiltrating immune cells and differential gene expression levels in various cancer types (Li et al., 2020). We used the TIMER database to investigate the differential expression of DLX1 in tumors and normal tissues in various cancer types.

RNA-sequencing data collection and bioinformatic analysis

We downloaded RNAseq data in level 3 HTSeq-FPKM format from 535 LUAD tumor tissues and 59 normal tissues from the TCGA (https://portal.gdc.cancer.gov/) LUAD project, as well as the corresponding clinicopathologic information. Data were transformed accordingly: fragments per kilobase per million (FPKM) RNAseq data were log2 transformed, and FPKM RNA sequencing data were converted to transcripts per million reads (TPM). The TCGA database is publicly available and written informed consent was obtained from patients prior to data collection. The 535 LUAD samples were divided into high and low DLX1 expression groups according to median DLX1 expression. DESeq2 was used with the absolute value of log-fold change > 1.5 and p-value < 0.05 as threshold parameters (Love, Huber & Anders, 2014). To analyze the DEGs between the two groups, we used the “ggplot2” software package to display the volcano map and heat map.

Functional enrichment analysis of DLX1-related DEGs in LUAD

Functional annotation and gene set enrichment analysis (GSEA) of DEGs were performed using the ClusterProfiler software package (Yu et al., 2012). Using absolute log fold change >1.5 and p < 0.05 as threshold parameters, we found 223 DEGs. The GSEA software was run using the MSigDB c2 collection (c2.all.v7.2.symbols.gmt). The functional enrichment analysis of DLX1-related DEGs in LUAD was visualized.

Kaplan–Meier plotter pan-cancer RNA-seq database analysis

The Kaplan–Meier mapper can be used to assess the prognostic impact of different genes in different types of cancer (http://kmplot.com). We performed the survival analysis of DLX1 expression in LUAD, as well as the prognostic value of DLX1 according to immune cells.

Correlation analysis of DLX1 expression level and immune cell infiltration in LUAD

The ssGSEA algorithm in GSVA (v1.34.0) R package (Hänzelmann, Castelo & Guinney, 2013) was used to evaluate the tumor penetration status of 24 kinds of immune cells (Bindea et al., 2013), including activated dendritic cells (aDC), B cells, CD8 T cells, cytotoxic cells, dendritic cells (DC), eosinophils, immature DC (iDC), macrophages, mast cells, neutrophils, CD56^bright natural killer (NK) cells, CD56dim NK cells, NK cells, plasmacytoid DC (pDC), T cells, T helper cells, T central memory (Tcm), T effector memory (Tem) cells, T follicular helper (TFH) cells, T gamma delta (Tgd) cells, Th1 cells, Th17 cells, Th2 cells, and Treg cells. The relationship between DLX1 expression levels and immune cell infiltration status was determined using Spearman correlation analysis.

DLX1 in CpG island DNA methylation status analysis

The methylation status of the DLX1 gene CpG locus in the TCGA-LUAD dataset was analyzed using MethSurv (https://biit.cs.ut.ee/methsurv/). Additionally, the relationship between DLX1 gene CpG methylation status and LUAD overall survival (OS) was evaluated.

Genetic changes in LUAD samples

The following four LUAD datasets were analyzed in cBioPortal for genomic alterations in the DLX1 gene (https://www.cbioportal.org/): Broad, cell 2012; MSK, JThorac Oncol 2020; OncoSG, Nat Genet 2020; and TCGA, Firehose Legacy. K-M survival curve analysis and log-rank test were used to determine the prognostic significance of DLX1 mutations.

Correlation analysis between expression level of DLX1 and clinicopathologic features in LUAD patients

We extracted clinicopathologic data from the TCGA-LUAD project and previously published studies in LUAD patients, including overall survival (OS), disease-specific survival (DSS), and progression-free interval (PFI) (Liu et al., 2018). Software analysis was used to compare the differences in various clinicopathologic parameters between the high and low DLX1 expression groups, such as T stage, N staging, M staging, DSS events, OS events, PFI events, and pathologic staging. We used a normality test (p < 0.05) for normally distributed data, and Kruskal–Wallis test and Dunn’s multiple methods to analyze intergroup differences. The Bonferroni method was used to correct the significance level results. The “ggplot2” software package was used for statistical visualization. Logistic regression analysis was used to evaluate the relationship between DLX1 expression levels and the clinicopathologic characteristics of LUAD patients.

Prognostic significance of DLX1 expression in LUAD

Survival data of LUAD patients from the TCGA-LUAD project and previously published data (Liu et al., 2018) were analyzed using the “Survival” package (statistical analysis of survival data) and the “Survminer” package (visualization). Prognostic analyses were performed. To determine LUAD survival according to DLX1 expression level, we performed K-M survival curve analysis and univariate and multivariate COX regression analyses. Diagnostic ROC curves and nomogram model analysis were performed using “proc,” “time ROC” (statistical analysis), and “ggplot2” (visualization) software packages to evaluate the predictive value of DLX1 expression level in LUAD diagnosis. The K-M survival curve was used for prognostic analysis of a subset of LUAD patients. Results included pooled sample sizes (percent), hazard ratios (HRs), confidence intervals (CIs), and P values. Forest plots were generated using the “ggplot2” software package.

Cell lines, cell culture, and siRNA transfection

A549, HCC-827, H1299, and SPC-A-1 were obtained from the Shanghai Institutes for Biological Sciences, China, and the BEAS-2B lung epithelium cell line was obtained from the cell bank of the Kunming Animal Research Institute. All LUAD cell lines were grown in RPMI-1640 medium (Corning, Corning, NY, USA) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin and streptomycin. BEGM medium (CC-3170, Lonza, Basel, Switzerland) was used to grow the BEAS-2B lung epithelial cell line. All cell lines were grown at 37 °C in humid air. H1299 and A549 cells were seeded in six-well plates 24 h before transfection. When the cell growth area was approximately 70%, cells were transfected with siRNA targeting DLX1 and negative control siRNA using Lipofectamine RNAimax reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer’s protocol. The efficiency of the transfection was assessed by quantitative real-time RT-PCR (qRT-PCR) analysis.

Cell proliferation and cell migration assay

After cells were transfected with siRNA to target specific genes and negative control siRNA for 24 h, cells were harvested and counted in a cell counting chamber. Cells were then plated on 96-well plates at a density of 2,000 cells per 100 ml well. Cell proliferation was detected using the Cell Counting Kit-8. The optical density (OD) was measured at 450 nm every 24 h and for four consecutive days after the addition of CCK8 solution according to the manufacturer’s instructions. LUAD cells transfected with si-DLX1 and si-NC were cultured in six-well plates for 24 h. The surface of the cell layer was scratched with a pipette tip. Images were taken with a microscope 0 h and 24 h after damage. The distance of the damaged area at 24 h was measured. The relative mobility was calculated and evaluated by comparing the distance of the damaged area at 0 h. Cell migration and invasion were detected using Transwell chambers (Corning). To perform cell migration assays, 2 × 10⁵ cells were added to the upper chamber and resuspended with 100 µL of serum-free RPMI 1640 medium, and 600 µL of RPMI 1640 medium containing 10% FBS was added to the lower chamber. After 24 h, the cells on the bottom surface of the lower chamber membrane were fixed with 4% paraformaldehyde and stained with 0.5% crystal violet, and the staining in five random fields was calculated by counting the number of cells under the microscope. For cell invasion, Matrigel was added to the upper chamber 3 h in advance. The remaining steps were performed as described above.

Real-time RT-PCR assay

The cells were lysed by RNA iso Plus (Takara Bio, Shiga, Japan) using real-time RT-PCR. This research use primer was as follows:

β-actin-F: GACCTGACTGACTACCTCATGAAGATCGC,

β-actin-R: GTCACACTTCATGATGGAGTTGAAGG,

DLX1-F: CGCTTCAATGGCAAGGGAAA,

DLX1-R: CTCTCCGGCAGAGCTAGGTA.

Statistical analysis

DLX1 expression in the normal and LUAD groups was statistically analyzed using the Wilcoxon rank sum test. We divided the patients into two groups according to the median expression of DLX1. Clinicopathological characteristics of DLX1 were analyzed using the Wilcoxon rank sum test, Kruskal–Wallis test, and logistic regression. Prognostic analysis was performed using Kaplan–Meier survival analysis and Cox univariate and multifactorial analysis. To evaluate the diagnostic significance of DEGs, ROC curves were plotted using the “plotROC” program package. The “RMS” R program package was used to plot the prognostic value of LUAD patients on the nomogram. All statistical analyses were performed in the R environment (v3.6.3). Wound healing and Transwell experiments were analyzed using Image J software. We used GraphPad Prism 8.0 software to graph the data, Student’s t-test to determine the significance of data between two experimental groups, and one-way ANOVA for multiple group comparisons. Each experiment was repeated at least three times and expressed as mean ± standard deviation (mean ± SD). P values < 0.05 were considered statistically significant differences (*: p < 0.05, **: p < 0.01, ***: p < 0.001).

DLX1 expression levels are significantly elevated in multiple cancers including LUAD

When analyzing the expression of DLX1 in different tumor types online in TIMER (https://cistrome.shinyapps.io/timer/), we found that the expression of DLX1 was significantly higher in malignant tissues than in normal tissues (p < 0.05, Fig. 1A). DLX1 expression in different malignant tumor types was analyzed using the TCGA database, followed by paired difference analysis. The results showed that DLX1 expression was significantly higher in tumor tissues compared to normal tissues (p < 0.05, Fig. 1C). We also determined that DLX1 expression in pan-cancer paired cancer tissues and adjacent normal tissues by means of the TCGA dataset. DLX1 levels were significantly higher in 12 of the 18 cancers compared to normal tissues (p < 0.05, Fig. 1B). DLX1 was significantly more expressed in unpaired samples (Fig. 1C) versus paired samples (Fig. 1B) in LUAD than in normal tissues (p <0.05).

The expression of DLX1 in LUAD was up-regulated, which indicated a poor prognosis in patients with LUAD

DLX1 expression and clinical significance were analyzed in the TCGA database. We found that OS, DSS, and PFI were shorter in LUAD patients with higher DLX1 levels (Figs. 2A–2C). We further investigated the diagnostic value of DLX1 in the differentiation of LUAD samples from normal lung tissue. The ROC curve analysis confirmed that the AUC value of DLX1 was 0.754, 95% CI [0.703–0.805] (Fig. 2D). These results suggest that DLX1 is upregulated in LUAD and that high expression levels of DLX1 are associated with poor prognosis in lung cancer patients. Based on the validation of the prognostic value of DLX1 in different subgroups, we further determined the prognostic value of DLX1 in different clinical subgroups, namely pathological stage, tumor lymph node metastasis (TNM) stage, gender, age, race, smoking status, and smoking age. The results showed that the upregulation of DLX1 levels was associated with poor clinical outcomes in lung cancer patients (Figs. 2E–2G).

Relationship between the expression of DLX1 and clinical characteristics of LUAD

We found that DLX1 expression was significantly correlated with pathological stage, OS events, DSS events, and PFI events in LUAD patients (Table 1). In addition, logistic analysis showed that DLX1 upregulation correlated with T stage (T2 & T3 & T4 vs. T1) and pathological stage (stage II, stage III, and stage IV vs. stage I) (Table 2).

Functional enrichment analysis of DLX1-related DEGs in LUAD

Based on the median DLX1 expression, the 535 LUAD patients were divided into high DLX1 expression group and low DLX1 expression group. We then found 223 DEGs in the high DLX1 expression group compared to the low DLX1 expression group (184 up-regulated and 39 down-regulated) using absolute log fold change >1.5 and p < 0.05 as threshold parameters (Fig. 3A). The single gene co-expression heat map in Fig. 3B shows the top 100 DEGs (Fig. 3B).We then performed functional annotation of DLX1-associated DEGs in LUAD patients using the “clusterProfiler” R package. The results of GO and KEGG enrichment analysis included highly enriched biological process (BP), cellular competent (CC), and molecular function (MF) (p < 0.05). The most common biological processes included “multicellular biological processes,” “apoptotic processes involved in development,” and “apoptotic processes involved in development” (p < 0.05). The most common biological processes included “multicellular biological processes,” “apoptotic processes involved in development,” “positive regulation of cell secretion,” and “γ-aminobutyric acid signaling pathway.” The most abundant cellular components were “GABA receptor complex,” “synaptic membrane,” “ion channel transporter complex,“ and “transport vesicles.” The most abundant cellular components were “GABA receptor complex,” “synaptic membrane,” “ion channel transporter complex,” “transport vesicles,” “cell base, and “intermediate filaments”. The most abundant molecular functions were “channel activity,” “motor activity,” “neurotransmitter receptor activity,” and “DNA-binding transcriptional repressor activity.” The most abundant molecular functions were “channel activity,” “motor activity,” “neurotransmitter receptor activity,” and “DNA binding transcriptional repressor activity.” The most common metabolic pathways include “ascorbate and aldehyde metabolism,” “chemical carcinogenesis,” “steroid hormone biosynthesis,” “neuroactive ligand–receptor interaction,” and “neurotransmitter receptor activity” (Figs. 3C–3F). GSEA revealed that the DEGs associated with DLX1 were significantly enriched in a cluster of cell proliferation-related genes (Fig. 4), including genes associated with cell cycle checkpoints (Nes = 2.012, Padj = 0. 036, FDR = 0.029), G2-M checkpoints (Nes = 1.780, Padj = 0.036, FDR = 0.029), M phase (Nes = 1.952, Padj = 0.036, FDR = 0.029), S phase (Nes = 1.589, Padj = 0.047, FDR = 0.037), mitotic G2-M phase (Nes = 1.867, Padj = 0.036, FDR = 0. 029), DNA repair (Nes =1.725, Padj = 0.036, FDR = 0.029), DNA replication (Nes = 1.752, Padj = 0.036, FDR = 0.029), gene expression regulation (Nes = 1.623, Padj = 0.042, FDR = 0.034), and DNA strand extension (Nes = 1.899, Padj = 0. 036, FDR = 0.029.) DLX1-associated DEGS were significantly associated with Fceri-mediated MAPK activation regulated by TP53 activity and MET activation of PTK2-regulated signaling pathways (Table S2).

The expression level of DLX1 is related to the infiltration of various immune cells in LUAD

ssGSEA was used to evaluate the infiltration of 24 types of immune cells in LUAD tissues. To evaluate the correlation between immune cell infiltration and expression, the association with DLX1 was evaluated using Spearman correlation analysis (Fig. 5A).The expression level of DLXI was positively correlated with Th2 cells (R = 0.167) (Fig. 6B) and negatively correlated with B cells (R = 0.091) (Fig. 5B), CD8 T cells (R = −0.137) (Fig. 5C), cytotoxic cells (R = −0.118)) (Fig. 5D), eosinophils (R = −0.149) (Fig. 5E), iDC (R = −0.089) (Fig. 5F), DC (R = −0.086) (Fig. 5G), Th17 cells (R = −0.132) (Fig. 6C), T cells (R = −0.110) (Fig. 6D), TFH (R =−0.094) (Fig. 6E), Mast cells (R = 0.157) (Fig. 6F), and neutrophils (R = −0.133) (Fig. 6G). The level of immune cell infiltration in the tumor was consistent with the results of Spearman analysis (Fig. 6A). Analysis of OS in immune cell infiltration DLX1 expression by Kaplan–Meier mapping (http://kmplot.com) showed poor DLX1 survival in high expressing type 2 helper T cells and poor DLX1 survival in low expressing B cells, CD4T cells, and eosinophils (Fig. 7).

Methylation status of DLX1 gene is associated with prognosis of patients with LUAD

DNA methylation levels in the DLX1 gene and the prognostic value of CpG islands in DLX1 gene were analyzed using the MethSurv tool. The results showed 32 methylated CpG islands, including cg07936950, cg05938001, cg12308746, and cg17737681, with elevated DNA methylation levels (Fig. 8). Furthermore, the methylation levels of two CpG islands (cg05938001 and cg06996609) were prognostic (p < 0.05) (Table 3). Higher DLX1 methylation in these two CpG islands, particularly cg05938001, was associated with worse OS in LUAD patients compared with LUAD patients with lower DLX1 CpG methylation.

The gene change of DLX1 had no effect on the survival outcome of patients with LUAD

Genetic alterations in the DLX1 gene were then analyzed using samples from 1,678 LUAD patients from the following four datasets: Broad, cell 2012 (n = 183); MSK, J Thorac Oncol 2020 (n = 604); OncoSG, Nat Genet 2020 (n = 305); and TCGA, Firehose Legacy (n = 586). Genetic alterations in the DLX1 gene were observed in only 1.5% of LUAD patients (Fig. 9A). K-M survival curves showed no significant difference in OS (p = 0.209) and DSS (p = 0.444) in patients with or without DLX1 gene alterations (Figs. 9B, 9C).

Univariate and multivariate Cox regression analysis of OS, DSS, and PFI with different parameters

We performed univariate Cox regression analysis in the TCGA-LAUD cohort to determine whether DLX1 expression level could be used as a valuable prognostic biomarker (Table 4). Univariable Cox regression results showed that high DLX1 expression, pathological stage, TNM stage, presence of residual lesions after treatment, and treatment outcome were associated with OS and DSS in LUAD patients (Figs. 10A–10B). To determine whether DLX1 expression level could be an independent prognostic factor for LUAD patients, multifactorial Cox regression analysis was performed. We confirmed that increased DLX1 expression was a significant independent prognostic factor in the TCGA-LUAD cohort, which was directly correlated with pathological stage, the presence of residual foci after treatment, and the outcome of treatment (Figs. 10A–10B). Cox regression analysis revealed that DLX1 expression, pathologic staging, T-stage, N-stage, presence of residual disease after therapy, and outcome were associated with PFI in patients with LUAD (Fig. 10C). To determine whether DLX1 expression level could be an independent prognostic factor for LUAD patients, multifactorial Cox regression analysis was performed. We found that DLX1 expression was a significant independent prognostic factor in the TCGA-LUAD cohort, with a direct correlation with pathologic stage, presence of post-treatment residual lesions, and response rate (Fig. 10C).

Construction and verification of normograph based on DLX1

DLX1 was found to be an independent prognostic factor for LUAD in a multifactorial analysis. By combining DLX1 expression with TNM stage, pathological stage, and presence of residual foci after treatment, we constructed predictive models for OS, DSS, and PFS. To incorporate DLX1 as a LUAD biomarker, we constructed a nomogram. Higher overall OS, PFS, and DSS scores were associated with poorer prognosis (Figs. 11A–11C).

Depletion of DLX1 significantly suppressed proliferation and migration of LUAD cells

To detect the expression of DLX1, we examined the expression level of DLX1 in LUAD cell lines using qRT-PCR. The results confirmed that DLX1 was significantly up-regulated in lung cancer cell lines, specifically H1299 and A549 cells (Fig. 12A). qRT-PCR analysis showed that DLX1 mRNA expression was significantly decreased after transfection of H1299 and A549 cells with si-DLX1 (Figs. 12B–12C). Cell counting kit-8 (CCK8), Transwell, and wound healing assay analysis showed that the deletion of DLX1 significantly inhibited the cell proliferation and the cell migration ability of LUAD (Figs. 12D–12K). These data indicate that DLX1 is highly expressed in LUAD and significantly affects its proliferation, migration, and invasion.