The expression of nine HOX genes of the HOXD cluster in LUAD tissues and normal tissues were examined using The University of Alabama at Birmingham Cancer data analysis Portal (UALCAN) database (http://ualcan.path.uab.edu/analysis-prot.html) (26). The prognostic value of HOXD1 expression for overall survival (OS), first progression (FP) and post-progression survival (PPS) was analyzed using the Kaplan-Meier plotter (http://kmplot.com) (27).

The human normal lung epithelial cell line (BEAS-2B) and LUAD cell lines (A549, H1299, H1650 and H1975) were purchased from Cell Bank/Stem Cell Bank, Chinese Academy of Sciences. BEAS-2B and A549 were cultured in DMEM (Shanghai VivaCell Biosciences, Ltd.) supplemented with 1% penicillin/streptomycin (Beijing Solarbio Science & Technology Co., Ltd.) and 10% FBS (CellMax). H1299, H1650 and H1975 cells were cultured in RPMI 1640 (Shanghai VivaCell Biosciences, Ltd.) supplemented with 1% penicillin/streptomycin (Beijing Solarbio Science & Technology Co., Ltd.) and 10% FBS. All cells were cultured at 37°C and 5% CO₂.

The lentiviruses used in this study were produced by Shanghai GeneChem Co., Ltd. Lentiviruses expressing HOXD1 were produced using the lentiviral vector GV492 (Shanghai GeneChem Co., Ltd.). The virus transfection process was performed according to manufacturer's protocol. Briefly, A549 or H1299 cells were seeded at 3–5×10⁴ cells/ml in 6-well plates. Lentiviruses containing the negative control (LV-NC; empty lentiviral vector GV-492) or HOXD1 vector (LV-HOXD1) were transfected into A549 or H1299 cells (multiplicity of infection, 10). The volume of culture medium was 1 ml/well in 6-well plates and 1X HitransG P (Shanghai GeneChem Co., Ltd.) was added to each well and incubate for 14–26 h at 37°C and 5% CO₂, the medium was replaced with fresh complete culture medium. Puromycin (6 µg/ml) selection was then used for 2 days and 3 µg/ml puromycin was used for maintenance to produce stably transfected cells.

A total of 16 female BALB/cA-nu mice (weight, 18–20 g; age, 4–5-weeks; Beijing HFK Bioscience Co., Ltd.,) were housed under specific pathogen-free conditions. The mice had ad libitum access to food/water and were housed under a 12 h light/dark cycle. The temperature was maintained at 24±2°C and a relative humidity range of 50–60% was maintained. All procedures were approved by the Animal Use and Care Committee at Shandong Provincial Hospital Affiliated with Shandong First Medical University (approval no. 2021-622; Jinan, China). A549 cells transfected with either the LV-NC or LV-HOXD1 plasmid vector via lentiviruses were used in the tumor formation assay. The nude mice were divided into two groups (8 mice/group): LV-NC and LV-HOXD1. All mice were injected subcutaneously in the right hind limb with 3×10⁶ A549 LV-NC or LV-HOXD1 cells in 150 µl PBS. The experimental duration of the present study was 5 weeks, with the first week allocated for adaptive feeding. The experiment commenced at the beginning of the second week, with the subcutaneous injection of cells to establish a xenograft tumor animal model. Animal health was observed daily and the body weight of the mice and tumor diameters were measured weekly for 4 weeks. If any humane endpoints were reached, the animals were sacrificed. These included a tumor diameter >20 mm, weight loss >20% of body weight, the animal exhibited cachexia or wasting syndrome or the size of the solid tumor >10% of body weight. Notably, none of the mice succumbed to humane endpoints during the experimental process. The mice were euthanized by cervical dislocation. Tumors were then dissected and tumor weights were measured.

The siRNAs (siR) used in the present study were constructed by external companies (siR-DNMTs, Guangzhou RiboBio Co., Ltd.; siR-NC, Shanghai GenePharma Co., Ltd.). A549 cells were incubated with 10 µM siR-DNMTs (Guangzhou RiboBio Co., Ltd.) and siR-negative control (siR-NC) (Shanghai GenePharma Co., Ltd.) using Lipofectamine^® RNAiMAX Reagent (cat. no. 13778-150; Invitrogen; Thermo Fisher Scientific, Inc.) (Table SI). An equal concentration of siR-NC was added to A549 cells as a control. Cells were seeded in 6-well plates until they reached 40–50% confluence. To prepare the transfection mix according to the manufacturer's protocol, 9 µl Lipofectamine^® RNAiMAX was diluted in 150 µl Opti-MEM Medium (Gibco; Thermo Fisher Scientific, Inc.) and mixed. Additionally, 3 µl siRNA was diluted in 150 µl Opti-MEM Medium and mixed separately. The diluted siRNA was then mixed with the diluted Lipofectamine^® RNAiMAX in a 1:1 ratio to form a siRNA-lipid complex. This complex was incubated for 5 min at room temperature before being added to the cultured cells. Following transfection, the cells were incubated at 37°C. After 48 h of transfection, the medium was removed from the wells. Cells were washed 3 times with cold PBS and then 500 µl RNAiso Plus reagent (Takara Bio, Inc.) was added to thoroughly lyse cells using RNase-Free Pipette Tips.

Total RNA was extracted from BEAS-2B cells and the LUAD cell lines, A549 cells transfected with siR-NC or with siR-DNMTs, and A549 and H1299 cells transfected with lentiviruses containing the negative control (LV-NC) or HOXD1 (LV-HOXD1) vector using RNAiso Plus reagent (Takara Bio, Inc.) according to the manufacturer's instructions. The RNA pellet was recovered in RNase-free water and the RNA concentration and purity were measured using the NanoDrop 1000 (Thermo Fisher Scientific, Inc.). cDNA was synthesized using the Evo M-MLV RT Premix (Hunan Accurate Bio-Medical Technology Co., Ltd.) according to the manufacturer's protocol. The cDNA was used as the template for RT-qPCR using the Taq SYBR Green qPCR Premix (Jiangsu Best-Enzymes Biotechnology Co., Ltd.). The following thermocycling conditions were used for qPCR: Initial denaturation at 95°C for 10 sec; 40 cycles of 60°C for 10 sec and 72°C for 30 sec. The final results were analyzed as previously described (28) and the relative quantification of target genes were analyzed using the 2^−ΔΔCq method (29) after normalization to β-actin or GAPDH expression levels. The primers used are listed in Table SII.

Cells were lysed on ice using RIPA lysis buffer (New Cell & Molecular Biotech Co., Ltd.) containing Protease Inhibitor Cocktail (New Cell & Molecular Biotech Co., Ltd.). Protein concentrations in the samples were quantified using the NanoDrop 1000 (Thermo Fisher Scientific, Inc.). A total of 75 µg protein/lane were separated by SDS-PAGE electrophoresis using a 10–12.5% gel, followed by transfer to PVDF membranes (cat. no. IPVH00010; MilliporeSigma). After blocking in 5% skim milk for 2 h at room temperature, the membranes were incubated with primary antibodies against HOXD1 (1:500; cat. no. ab220856; Abcam), Flag M2 (1:5,000; cat. no. m20008m; Abmart Pharmaceutical Technology Co., Ltd.) and GAPDH (1:50.000; cat. no. 60004-I–Ig; Proteintech Group, Inc.) overnight at 4°C. After washing three times with TBST containing 0.1% Tween-20, the membranes were incubated with anti-mouse IgG HRP-linked antibodies (1:20,000; cat. no. AB0102; Shanghai Abways Biotechnology Co., Ltd.) and anti-rabbit IgG HRP-linked antibodies (1:50,000; cat. no. AB0101; Shanghai Abways Biotechnology Co., Ltd.) at room temperature for 1 h and an ECL kit (cat. no. WBKLS0500; MilliporeSigma) were used for visualization. The chemiluminescence imaging system (ChemiDoc M, Bio-Rad Laboratories, Inc.) was used for imaging and ImageJ software (version 1.51j8; National Institutes of Health,) was used for semi-quantitation of protein levels. GAPDH was used as a reference control.

The UALCAN database was used to analyze HOXD1 promoter methylation levels (26). CpG islands in the HOXD1 promoter sequence were predicted using MethPrimer (http://www.urogene.org/methprimer/) (30). TBS was performed by Igenebook Biotechnology Co., Ltd. Briefly, bisulfite sequencing PCR primers targeting the HOXD1 promoter region were designed using MethPrimer (Table SIII). The HOXD1 promoter region was selected from upstream 2,000 bp to downstream 1,000 bp of the transcription start site (TSS) and located on the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/) nucleotide sequence starting from 176,186,579 and ending at 176,189,579 (accession no. NC_000002). Genomic DNA from BEAS-2B, A549 and H1299 cells were extracted and bisulfite treatment was performed using the EZ DNA Methylation-Gold^™ Kit (cat. no. D5005; Zymo Research Corp.). DNA quality was detected using the Qubit 4.0. Bisulfite-treated templates were subjected to bisulfite sequencing PCR (BSP) amplification by high-fidelity U-base-resistant DNA polymerase BSP amplification. BSP amplification products from the same sample were mixed and labeled primers were amplified to obtain a bisulfite-converted DNA library. The library from each sample was pooled and sequenced using the Illumina NovaSeq6000 platform (Illumina, Inc.) using the paired end 150 bp method as previously described (31,32). A library concentration of 13.1 ng/µl was measured using Qubit 4.0. Trimmomatic (version 0.36), BSMAP (version v2.7.3) and FastQC (version 0.11.7; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) were used for data analysis (33,34).

A549 cells at a density of 70–80% were treated with 0.1, 1.0, 5.0 or 10.0 µM decitabine [5-aza-2′-deoxycytidine (DAC); cat. no. ID0120; Beijing Solarbio Science & Technology Co., Ltd.] and 10 µM DMSO (cat. no. D8371; Beijing Solarbio Science & Technology Co., Ltd.) for 96 h at 37°C and 5% CO₂. The cells were then collected for RNA extraction and RT-qPCR analysis.

Cell proliferation was measured using crystal violet staining, Cell Counting Kit-8 (CCK-8) and colony formation assays. For the crystal violet staining assay, cells were seeded in a 12-well culture plate at 1×10⁵ cells/well and incubated for 3 days. After washing twice with cold PBS, cells were fixed using 5% glacial acetic acid for 15 min at room temperature and then stained with 0.1% crystal violet for 10 min at room temperature. Cells were imaged using a camera. For the CCK-8 assay, it was performed according to manufacturer's protocol of CCK-8 reagent (cat. no. BS350B; Biosharp Life Sciences). Briefly, cells were seeded in a 96-well culture plate at a density of 2–3×10³ cells/well. After the cells were adherent and were cultured for different time periods (0, 24, 48 and 72 h), 10 µl/well of CCK-8 reagent was added, and then the cells were incubated for 2 h at 37°C and 5% CO₂. Absorbance at 450 nm was measured using a microplate reader (SpectraMax 190; Molecular Devices, LLC). For the colony formation assay, cells were seeded at a density of 1×10³ cells/well in a 6-well culture plate and incubated for 7 days until a single cell proliferated to form a visible cluster which was defined as a colony. After discarding the culture medium and washing twice with PBS, 4% paraformaldehyde was added to each well to fix the cells for 30 min at room temperature. Then, cells were stained with 2 ml of 0.1% crystal violet for 3 min at room temperature. Visible colonies were counted using Image J software.

For the wound-healing assay, 1×10⁶ A549 cells were seeded in a 6-well plate and then scratched with a 10 µl pipette tip. The cells were washed 3 times with PBS to remove the scratched cells and then serum-free medium (DMEM or RPMI 1640) was used to culture the cells. Cells were cultured at 37°C and 5% CO₂. The migration of cells at the indicated time points of 0 and 96 h was observed using an inverted fluorescence microscope (Olympus IX73; Olympus Corporation). The width of the wound was measured using Photoshop (version 2017.1.6; Adobe Systems, Inc.). For the Transwell experiments, Transwell inserts were used with or without a Matrigel coating (cat. no. HY-K6001; MedChemExpress) for the cell invasion and migration experiments. Matrigel was melted at 4°C overnight and diluted with pre-cooled serum-free medium at 4°C to a final concentration of 1 mg/ml and maintained on ice. Then, 100 µl of diluted Matrigel was added to the center of the bottom of the upper chamber and incubated at 37°C for 4–5 h to dry. Control cells (A549 LV-NC and H1299 LV-NC) and HOXD1-overexpression cells (A549 LV-HOXD1 and H1299 LV-HOXD1) were seeded at a density of 2×10⁵ cells in the upper chamber in serum-free medium. Regular culture medium containing 10% FBS was added to the lower chamber. Cells were then incubated at 37°C for 24 h. Migratory and invasive cells in the lower chambers were stained using 0.1% crystal violet for 2 min at room temperature. Finally, cells were imaged using a light microscope (Olympus IX73; Olympus Corporation). Cells were counted and photographed in three randomly selected fields of view and quantified using Image J software.

The ChIP assays were conducted by Igenebook Biotechnology Co., Ltd. Briefly, ~3×10⁷ A549 cells were washed twice in cold PBS, cross-linked with 1% formaldehyde for 10 min at room temperature, quenched with glycine for 5 min at room temperature and then washed twice with cold PBS at room temperature. Cells were collected by centrifugation at 1,000 × g for 5 min at 4°C. Samples were lysed using 50 mM Tris-HCl (pH 8.0), 10 mM EDTA, 1% SDS, 1X protease inhibitor cocktail (cat. no. 5056489001; MilliporeSigma) and chromatin on ice. The chromatin was sheared into an average DNA fragment length of 200–500 bp. Additionally, 20 µl of chromatin was stored at −20°C for input DNA and 100 µl chromatin was incubated at 4°C overnight with antibodies against Flag M2 (1:50; cat. no. 14793; Cell Signaling Technology), DNMT1 (1:100; cat. no. 24206-1-AP; Proteintech), DNMT3A (1:100; cat. no. 20954-1-AP; Proteintech), DNMT3B (1:50; cat. no. 26971-1-AP; Proteintech) or IgG (1:100; cat. no. 2729S; Cell Signaling Technology, Inc.) as a negative control for immunoprecipitation. Then, 30 µl of protein beads (Dynabeads^™ protein G; cat. no. 10004D; Thermo Fisher Scientific, Inc.) were added and the samples were further incubated for 3 h at 4°C. The beads were then washed once with 20 mM Tris/HCL (pH 8.1), 50 mM NaCl, 2 mM EDTA, 1% Triton X-100 and 0.1% SDS, washed twice with 10 mM Tris/HCL (pH 8.1), 250 mM LiCl, 1 mM EDTA, 1% NP-40 and 1% deoxycholic acid and twice with 1X TE buffer (10 mM Tris-Cl at pH 7.5 and 1 mM EDTA). Bound material was then eluted from the beads using 300 µl of elution buffer (100 mM NaHCO3 and 1% SDS), treated with RNase A (8 µg/ml; cat. no. EN0531; Thermo Fisher Scientific, Inc.) for 6 h at 65°C and then with proteinase K (345 µg/ml; cat. no. P6556; MilliporeSigma) overnight at 45°C. DNA concentration and purity were detected using the Q-bit (Qbit 4.0; Invitrogen; Thermo Fisher Scientific, Inc.). The purified products were used to construct sequencing libraries following the protocol provided by the I NEXTFLEX^® ChIP-Seq Library Prep Kit for Illumina^® Sequencing (cat. no. NOVA-5143-02, Bioo Scientific). The concentration of the libraries were assayed using the Qubit 4.0 and the fragment size determined using the QSep400 (Bioptic). The library concentration of 22.5 nM was determined by qPCR. Then, the purified products were sequenced on an Illumina NovaSeq using the paired end 150 bp method (35). Trimmomatic (version 0.36) was used to filter out low-quality reads. Next, clean reads were mapped to the human genome using BWA (version 0.7.15-r1140). Samtools (version 1.3.1) was used to remove potential PCR duplicates. MACS2 software (version 2.1.1.20160309) was adopted to screen peaks (bandwidth 300 bp; model fold 5, 50; q value 0.05). Peaks were assigned to genes if their midpoint was closest to the TSS of a single gene (36). HOMER (version 3) was used to predict motif occurrences within peaks with default settings for a maximum motif length of 12 bp (37). The ClusterProfiler (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) R package (38) was employed to perform Gene Ontology (GO; http://geneontology.org/) (39) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) (40) enrichment analyses. The GO and KEGG enrichment analyses were calculated using hypergeometric distribution with a q value cutoff of 0.05.

For ChIP-qPCR, the extracted DNA fragments were used as templates for qPCR analysis using the ChamQ SYBR Color qPCR Master Mix (cat. no. Q411; Vazyme Biotech Co., Ltd.). The 2,000 bp region upstream of the TSS in the HOXD1 promoter region was divided into seven segments (F1, F2, F3, F4, F5, F6 and F7), and primers were designed for each segment. Primers are listed in Table SIV. The following thermocycling conditions were used for qPCR: Initial denaturation at 95°C for 30 sec; 40 cycles of 95°C for 10 sec and 60°C for 30 sec. The data were normalized to the input. The final results were analyzed as previously described (28) and the relative quantification of target genes were analyzed using the 2^−ΔΔCq method (29).

Differential expression data for the HOXD family in lung adenocarcinoma derived from the UALCAN database were analyzed using Welch's t-test. Survival analyses from the Kaplan-Meier Plotter database were conducted using Cox regression analysis. All experiments were repeated at least three times. GraphPad Prism software (version 8.0.1; Dotmatics) was used for all statistical analyses. A student's t-test was used to analyze the colony formation, Transwell and xenograft tumor assays and the differential expression of BMP2 and BMP6. The CCK-8 assay data were analyzed using a two-way ANOVA followed by Bonferroni's post hoc test. HOXD1 expression in lung adenocarcinoma cell lines and data from the DAC treatment experiments were analyzed using one-way ANOVA followed by Dunnett's post hoc test. The ChIP-qPCR results were analyzed using a two-tailed unpaired Student's t-test. Data were presented as mean ± SD. P<0.05 was considered to indicate a statistically significant difference.

The UALCAN database was used to examine the expression levels of nine HOX genes from the HOXD cluster in LUAD (Fig. 1). It was demonstrated that HOXD1 expression level in the tumor group (median level, 3.125) was significantly lower compared with that in the normal group (median level, 3.735) (P<0.05), while HOXD3 (median level of tumor group, 0.211; median level of normal group, 0.084) and HOXD4 (median level of tumor group, 0.317; median level of normal group, 0.097) expression levels were significantly increased in the tumor groups compared with the normal groups (P<0.05). There were no significant differences in HOXD8 (median level of tumor group, 0.903; median level of normal group, 0.864) (P>0.05), HOXD9 (median level of tumor group, 0.582; median level of normal group, 0.47) (P>0.05) and HOXD10 (median level of tumor group, 0.082; median level of normal group, 0.049) (P>0.05) expression levels between the tumor and normal groups. There were no data for the expression of HOXD11, HOXD12 and HOXD13 in the samples analyzed. In addition, HOXD1 expression levels were demonstrated to be low in LUAD cell lines compared with a normal lung epithelial cell line (Fig. 2A and B).

The prognostic values of the HOXD1 gene in patients with LUAD were analyzed using Kaplan-Meier plotter. It was demonstrated that patients with LUAD who had lower expression levels of HOXD1 had a poorer predicted OS (P=0.0014) and FP (P=0.00023) compared with patients with high expression levels of HOXD1, whereas there was no difference in PPS between the two patient groups (P=0.4) (Fig. 2C). These data indicated that HOXD1 expression was downregulated in LUAD in comparison to the normal lung, and low HOXD1 expression levels may be associated with a poor prognosis in patients with LUAD.

To explore the function of HOXD1 in LUAD, two stably transfected cell lines (A549 and H1299) overexpressing HOXD1 were established. A549 and H1299 cells were transfected with LV-NC and LV-HOXD1 to establish HOXD1 overexpression cell lines (Fig. 3A). Crystal violet staining, CCK-8 assays and colony formation assays demonstrated that the upregulation of HOXD1 significantly inhibited cell proliferation compared with control cells (P<0.05; Fig. 3B-D). Additionally, to determine cell motility, a wound healing assay was performed. These results demonstrated that HOXD1 overexpression in LUAD cells was associated with a slower wound closure (Fig. 3E). Moreover, Transwell assays showed that the migration and invasion of LUAD cells were inhibited when HOXD1 was upregulated compared with control cells (Fig. 3F).

A549 LV-NC and A549 LV-HOXD1 cells were subcutaneously injected into nude mice and after 28 days, the mice were sacrificed and tumors were dissected. It was demonstrated that the weight of subcutaneous tumors in the HOXD1-upregulated group was lower compared with that in the control group (Fig. 3G). These assays indicated that HOXD1 may potentially have the ability to suppress the progression of LUAD.

To investigate the molecular mechanism of low HOXD1 expression levels in LUAD cells, the UALCAN database was used to analyze the methylation level of the HOXD1 promoter region in LUAD samples and normal samples. These results showed that the methylation level of LUAD samples was significantly higher compared with that of normal samples (Fig. 4A). This indicated that the expression of HOXD1 may be regulated by DNA methylation. CpG islands in the promoter region of the HOXD1 gene were predicted using MethPrimer. These results showed three CpG islands in the upstream 2,000 bp to downstream 1,000 bp of the TSS (Fig. 4B). Furthermore, the methylation level of the HOXD1 promoter in BEAS-2B, A549 and H1299 cells was measured, and the TBS results showed that the total level of methylation was not significantly different in A549 and H1299 cells compared with BEAS-2B cells (Fig. 4C). It is noteworthy that one CpG island in the posterior segment was significantly less methylated, and two CpG islands in the anterior section of this promoter region were significantly more methylated in LUAD cells compared with BEAS-2B cells (Fig. 4D). These findings suggested that in LUAD, the HOXD1 promoter region was regulated by DNA methylation.

The regulation mechanism of HOXD1 hypermethylation was investigated as it was hypothesized that DNMTs were involved in regulating HOXD1 expression in LUAD cells. To determine whether DNMTs were upstream regulators of HOXD1 expression in LUAD cells, A549 cells were treated with the DNMT inhibitor decitabine for 96 h. The mRNA expression level of the HOXD1 gene was detected using RT-qPCR, which demonstrated that DAC significantly increased HOXD1 expression levels in a dose-dependent manner (Fig. 5A). In addition, A549 cells were transfected with siRNAs targeting DNMT1, DNMT3A or DNMT3B and a significant increase in HOXD1 expression levels were demonstrated compared with negative controls (Fig. 5B). Furthermore, the 2,000 bp region upstream of the TSS in the HOXD1 promoter region was divided into seven segments, and primers were designed for each segment (Fig. 5C). ChIP-qPCR assays using A549 cells showed that DNMTs bound to the HOXD1 promoter region. DNMT1 combined with F1, F5 and F7 segments of the HOXD1 promoter (Fig. 5D). DNMT3A combined with F1, F2, F3, F4, F5, F6 and F7 segments of the HOXD1 promoter (Fig. 5E). DNMT3B combined with F1, F2, F3, F4, F5, F6 and F7 segments of the HOXD1 promoter (Fig. 5F). These ChIP-qPCR analyses suggested that the HOXD1 promoter region may be enriched with DNMT protein binding. These results suggested that DNMTs target the HOXD1 promoter region and could potentially contribute to local hypermethylation, inducing the dysregulation of HOXD1 in LUAD.

To better identify the downstream target gene potentially regulated by HOXD1 as a transcription factor in LUAD, ChIP-seq was performed on HOXD1 overexpressing A549 cells to analyze the downstream regulatory network of HOXD1. GO annotation analysis demonstrated that the target sequences of HOXD1 were most enriched in biological processes (Fig. 6A). After the GO enrichment data were sorted according to the P-value, from low to high, the top 20 GO-enriched functions were presented as a bubble graph, which reflected the number of genes and the degree of enrichment of HOXD1-enriched annotations to GO functions. (Fig. 6B). Moreover, following a descending order of P-value for all KEGG enrichment data, bar graphs representing the top 20 KEGG-enriched signaling pathways were created to show the enrichment of HOXD1 in various metabolic pathways (Fig. 6C).

Based on the results of ChIP-seq analysis, 24 genes associated with cellular processes were screened. Then, the expression of these genes was confirmed using RT-qPCR. These results showed that overexpression of HOXD1 elevated the mRNA expression levels of BMP2 and BMP6 (Fig. 6D).