DNA methylation-mediated differential expression of DLX4 isoforms has opposing roles in leukemogenesis

Previously, we reported the expression of DLX4 isoforms (BP1 and DLX7) in myeloid leukemia, but the functional role of DLX4 isoforms remains poorly understood. In the work described herein, we further determined the underlying role of DLX4 isoforms in chronic myeloid leukemia (CML) leukemogenesis. The expression and methylation of DLX4 isoforms were detected by real-time quantitative PCR (RT-qPCR) and real-time quantitative methylation-specific PCR (RT-qMSP) in patients with CML. The functional role of DLX4 isoforms was determined in vitro and in vivo. The molecular mechanism of DLX4 isoforms in leukemogenesis was identified based on chromatin immunoprecipitation with high-throughput sequencing (ChIP-Seq)/assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-Seq) and RNA sequencing (RNA-Seq). BP1 expression was increased in patients with CML with unmethylated promoter, but DLX7 expression was decreased with hypermethylated promoter. Functionally, overexpression of BP1 increased the proliferation rate of K562 cells with S/G2 promotion, whereas DLX7 overexpression reduced the proliferation rate of K562 cells with G1 arrest. Moreover, K562 cells with BP1 overexpression increased the tumorigenicity in NCG mice, whereas K562 cells with DLX7 overexpression decreased the tumorigenicity. Mechanistically, a total of 91 genes including 79 messenger RNAs (mRNAs) and 12 long noncoding RNAs (lncRNAs) were discovered by ChIP-Seq and RNA-Seq as direct downstream targets of BP1. Among the downstream genes, knockdown of RREB1 and SGMS1-AS1 partially revived the proliferation caused by BP1 overexpression in K562 cells. Similarly, using ATAC-Seq and RNA-Seq, a total of 282 genes including 151 mRNA and 131 lncRNAs were identified as direct downstream targets of DLX7. Knockdown of downstream genes PTPRB and NEAT1 partially revived the proliferation caused by DLX7 overexpression in K562 cells. Finally, we also identified and validated a SGMS1-AS1/miR-181d-5p/SRPK2 competing endogenous RNA (ceRNA) network caused by BP1 overexpression in K562 cells. The current findings reveal that DNA methylation-mediated differential expression of DLX4 isoforms BP1 and DLX7 plays opposite functions in leukemogenesis. BP1 plays an oncogenic role in leukemia development, whereas DLX7 acts as a tumor suppressor gene. These results suggest DLX4 as a therapeutic target for antileukemia therapy.

suppressor gene. These results suggest DLX4 as a therapeutic target for antileukemia therapy.

Background
Hematopoiesis is a process that is strictly controlled by gene expression regulation that drives differentiation from hematopoietic stem/progenitor cells to mature hematopoietic cells [1]. During leukemogenesis, biological processes including cell proliferation, expansion, self-renewal, and differentiation are disrupted, resulting in an accumulation of immature blast cells [2]. Hematological myeloid malignancies are a group of disorders characterized by clonal expansion of hematopoietic stem/progenitor cells, including acute myeloid leukemia (AML), chronic myeloid leukemia (CML), myelodysplastic syndromes (MDS) and myeloproliferative neoplasms (MPN). Epigenetic processes such as DNA methylation and histone modifications with roles in regulating gene(s) expression have been proved to play a vital role in myeloid malignancies [3]. Recurrent somatic mutations have been identified in genes involved in epigenetic regulators such as TET2, DNMT3A, IDH1/2 and EZH2 that were required for malignant transformation [4,5]. Furthermore, aberrant DNA hypermethylation of tumor suppressor genes is a frequent event and contributes to luekemogenesis in myeloid malignancies, and could also predict disease progression and clinical outcome [6]. Importantly, demethylation agents including azacytidine and decitabine have been approved as first-line therapy in AML and MDS [7].
The human DLX4 gene is a member of the DLX subfamily of homeobox genes, which has three RNA splicing isoforms encoding different proteins, namely DLX4 isoform 1 (NM_138281, encoding BP1), DLX4 isoform 2 (NM_001934, encoding DLX7), and DLX4 isoform 3 (encoding an unconfirmed protein). Most previous research on this topic has mainly focused on DLX4 isoform 1, thus BP1 is also called DLX4 in some studies. It has been demonstrated that DLX4 regulates diverse developmental processes such as skeletal patterning, neurogenesis, and hematopoiesis [8]. Moreover, BP1 is also implicated in biological progresses and early development, which are frequently dysregulated in human cancers, possible owing to DNA amplification [9,10]. BP1 is absent from most normal tissues but is commonly expressed in diverse human cancers including breast cancers and ovarian cancer, lung cancer, and acute leukemia [11][12][13][14][15][16][17]. Moreover, high BP1 expression promotes tumor progression and predicts clinical outcome in various human cancers [15][16][17][18]. These results suggest an oncogenic role of BP1 in human cancers, including leukemia.
Interestingly, epigenetic inactivation of DLX4 by DNA methylation has been proved in chronic lymphocytic leukemia (CLL) and several other solid tumors [19][20][21][22][23]. Moreover, DLX4 hypermethylation, with a role in silencing DLX4 expression, is associated with disease evolution in uterine cervical low-grade squamous intraepithelial lesions (LSILs) and non-small cell lung cancer (NSCLC) [22,23]. Importantly, our previous studies also confirmed the DLX4 hypermethylation pattern in myeloid malignancies [24][25][26], associated with poor prognosis in MDS and AML [24,25]. The results support the tumor suppressive role of DLX4 in human cancers, including myeloid malignancies. Superficially, there seems to be a contradiction regarding the role of DLX4 in tumorigenesis. Notably, we found that DLX4 hypermethylation was correlated with reduced DLX7 expression but not BP1 expression in AML and CML [25,26]. At the same time, our previous study also reported that BP1 expression is upregulated in AML, whereas DLX7 expression is downregulated [27]. These results suggest that BP1 and DLX7 may play opposite roles in luekemogenesis.
In the work described herein, we first evaluated the expression of two DLX4 isoforms, viz. BP1 and DLX7, in CML patients and investigated the potential correlations with clinical pathological parameters. Moreover, we determined the biological roles of BP1 and DLX7 in luekemogenesis by using in vivo and in vitro experiments. Finally, the potential molecular mechanism of BP1 and DLX7 in leukemogenesis was further explored. Collectively, the findings of the current study may provide broader insights into the understanding of leukemogenesis and provide potential therapeutic targets for antileukemia therapies.

Patients
The current investigation was approved by the Ethics Committee of the Affiliated People's Hospital of Jiangsu University (no. K-20190016, date 25/02/2019) in accordance with the Declaration of Helsinki, and included a total of 37 healthy controls and 75 CML patients who were treated at our hospital. Clinical and laboratory characteristics of CML patients are presented in Table 1. Bone marrow (BM) samples were obtained from all participants after written informed consent was signed. BM mononuclear cells were separated through density-gradient centrifugation by using lymphocyte separation medium (Solarbio, Beijing, China).

BSP
Bisulfite sequencing PCR (BSP) was further applied to investigate and/or confirm BP1 methylation by using TaKaRa Taq Hot Start Version (Tokyo, Japan) with the primers presented in Additional file 1: Table S1 [31]. Clone sequencing of PCR products was conducted as per previous literature [24][25][26]. Six independent clones from each sample were sequenced by Sanger sequencing (BGI, Shanghai, China).

Cell proliferation and cell cycle assays
Cell proliferation and cell cycle assays were carried out by cell counting and flow cytometry, respectively. For cell proliferation analysis, cells (1 × 10 5 cells/ml) were seeded into a six-well plate with complete medium for culture. After culturing for 0, 1, 2, and 3 days, cells were counted in a counting chamber by using an ordinary microscope (manual counting) for three times. For cell cycle analysis, cells (2 × 10 5 cells/ml) were seeded into a six-well plate with complete medium for culture. The BD Cycletest Plus DNA reagent kit (BD Pharmingen, San Diego, CA) was used to analyze the cell cycle distribution according to the manufacturer's protocols, followed by analysis using flow cytometry. Each experiment was repeated three times.

Xenograft mouse model
All mouse experiments were approved by the Committee on the Ethics of Animal Experiments of Jiangsu University (no. UJS-IACUC-AP-20190305073, date: 2019.03.05), also in compliance with the Basel Declaration. Mice of severe immune-deficient strain NCG (NOD/ShiLtJGpt-Prkdc em26Cd52 Il2rg em26Cd22 /Gpt) were obtained from GemPharmatech Co, Ltd (Nanjing, China). A total of 1 × 10 6 tested cells (K562-NC/K562-BP1/ K562-DLX7) were injected into each group of 5-week-old female NCG mice through the tail vein respectively. Body weight and peripheral blood of the mice were determined weekly. Growth of the mice affected by the leukemia cells was monitored (IVIS) at the 4th and 6th week. Mice were euthanized when they developed a bowed back and hind limb paralysis.

RNA-Seq
The details of RNA-Seq were reported in previous studies [32]. The sequencing data were analyzed with the assistance of Genesky Biotechnologies Inc. (Shanghai, China). Differentially expressed genes were analyzed by using Deseq2 software with P < 0.05 and |log 2 (fold change)| > 1.

ChIP-Seq
Chromatin immunoprecipitation (ChIP) assays were carried out by using the EZ-Zyme chromatin prep kit (Merck Millipore, Billerica, MA) based on the recommended protocols. Anti-BP1 antibody (NB100-481) (Novus Biologicals, Littleton, CO) and normal immunoglobulin G (IgG) were used for immunoprecipitation. Eluted DNA fragments were analyzed by high-throughput sequencing (ChIP-Seq) performed at Genesky Biotechnologies Inc. (Shanghai, China). All reads were paralleled to the National Center for Biotechnology Information (NCBI) human reference genome build 37 (GRCh37/ assembly hg19) using Bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml). The HOMER and MEME software packages were used to conduct a search of consensus binding motifs of BP1.

ATAC-Seq
Assay for transposase-accessible chromatin (ATAC) assay with high-throughput sequencing (ATAC-Seq) allows high-throughput sequencing of open chromatin regions with the help of transposases, which is similar to ChIP-Seq. The tested cells were processed to extract the nucleus, followed by transposition reaction. Separated and purified DNA was analyzed by ATAC-Seq performed at Genesky Biotechnologies Inc. (Shanghai, China). All reads were aligned to the NCBI human reference genome build 37 (GRCh37/ assembly hg19) using Bowtie2 (http:// bowtie-bio. sourc eforge. net/ bowti e2/ index. shtml). The HOMER and MEME software packages were applied to conduct a search of consensus binding motifs of DLX7.

FISH
Fluorescence in situ hybridization (FISH) assays were applied to determine the location of lncRNAs. Briefly, the exponential phase of cells seeded on slides were cleaned with phosphate-buffered saline (PBS) and fixed in 4% paraformaldehyde. Cells were then treated with the Cy3-labeled SGMS1-AS1 probe and incubated at 37 °C overnight. Moreover, Cy3-labeled U6 and 18S were also included as controls. The probe sequences are listed in Additional file 1: Table S1. After washing with hybridization solution and PBS for 5-10 min, cells were counterstained with 4′,6-diamidino-2-phenylindole (DAPI). Photos were captured by using an Olympus confocal laser scanning microscope.

Luciferase reporter assay
The putative binding area of human SGMS1-AS1 and the 3′-untranslated region (UTR) of human SRPK2 together with the matched mutant fragment were introduced individually into a GV268 vector (Genechem, Shanghai, China). Moreover, miR-181d-5p together with the matched mutant fragment was individually inserted into a GV272 vector (Genechem, Shanghai, China). To perform the reporter experiments, the constructed GV268 plasmid (mRNA-3′UTR or LncRNA-binding area) and constructed GV272 plasmid (miRNA) together with Renilla luciferase reporter plasmid were cotransfected into 293T cells by using X-tremeGENE HP transfection reagent (Roche, Basel, Switzerland). After 48 h of transfection, the luciferase activities and the firefly luciferase vitalities were examined under the Dual-Luciferase reporter assay system (Promega, Madison, WI). Relative firefly luciferase activity was normalized to Renilla luciferase activity to obtain the transfection efficiency.

RIP
RNA immunoprecipitation (RIP) assays were performed using a Magna RNA-binding protein immunoprecipitation kit (Millipore, Bedford, MA) on the basis of the recommended instructions. Briefly, 2 × 10 7 K562 cell lysates were incubated with RIP buffer containing magnetic beads conjugated with mouse IgG (negative controls) or human anti-Ago2 antibody (ab32381; Abcam, Cambridge, MA). The specimens were then incubated with proteinase K to isolate the immunoprecipitated RNA. RNAs were extracted and purified, then used for reverse transcription. Finally, RT-qPCR was carried out to confirm the presence of the binding targets. The primers given in Additional file 1: Table S1 were used to detect the SGMS1-AS1 and miR-181d-5p level.

Statistics
SPSS 22.0 and GraphPad Prism 5.0 were utilized for statistical analysis. Comparison of continuous variables between two groups was conducted by Independent T/Mann-Whitney U-tests, whereas comparison of categorical variables was performed by Pearson chi-squared/Fisher exact tests. The receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC) were applied to define the capability of BP1 and DLX7 expression for distinguishing CML patients from normal controls. In all analyses, (two-tailed) P values less than 0.05 were defined as statistically significant.

DNA methylation-mediated differential expression of DLX4 isoforms in CML
Following our previous study [27], we determined the expression pattern of DLX4 isoforms in CML patients. Interestingly, BP1 expression was significantly increased in CML patients (P = 0.006), whereas DLX7 expression was markedly decreased (P = 0.001) (Fig. 1a). To support these results, we further investigated the promoter methylation pattern of DLX4 isoforms in CML. Notably, independent CpG islands were identified in the promoter region of each DLX4 isoform (Fig. 1b). As reported in our previous study, the CpG island located at the promoter region of DLX7 (CpG island 2 in Fig. 1b) was methylated in CML patients and K562 cell line, and was associated with DLX7 expression [26]. However, the CpG island methylation located at the promoter region of BP1 (CpG island 1 in Fig. 1b) was undetectable in CML patients and K562 cell line by RT-qMSP, demonstrating that the CpG island located at the promoter region of BP1 was almost unmethylated in CML patients and K562 cell line. Moreover, the CpG island methylation located at the promoter region of BP1 was further validated by BSP; the results for K562 cell line and a representative CML patient are shown in Fig. 1c.
To investigate the clinical significance of DLX4 isoforms in CML, we first performed ROC curve analysis to evaluate the ability of DLX4 expression to discriminate AML from controls. The AUC value for BP1 and DLX7 were 0.624 (95% CI 0.520-0.727, P = 0.034, Fig. 1d) and 0.699 (95% CI 0.605-0.793, P = 0.001, Fig. 1e), suggesting that they may serve as potential biomarkers for distinguishing CML patients from controls. Based on ROC analysis, BP1 expression of 0.0256 and DLX7 expression of 0.0266 were set as cutoff points because the "sensitivity + specificity -1" reached the highest value. At the cutoff value, the sensitivity and specificity for BP1 expression were 33.3% and 97.3%, whereas they were 52% and 97.3% for DLX7 expression. According to the cutoff points, we divided the CML cases into two groups to analyze the clinical significance of both BP1 and DLX7 expression. No remarkable differences were found between groups with low and high BP1 expression among all clinical parameters (P > 0.05, Table 1). However, significant differences were observed between the groups with low and high DLX7 expression in the distribution of clinical . ∆CT reflects the disparity in CT value between control and target or reference sequences. The bone marrow sample from one normal control that possessed the minimal ∆CT between BP1/DLX7 and ABL1 transcript was selected as control and defined as 100% expression for BP1/DLX7 transcript. The median level of BP1/DLX7 expression in each group is shown by a horizontal line. b The coordinates of CpG islands in DLX4 gene. c Methylation density of BP1 promoter CpG island (CpG island 1) in a representative CML patient and K562 cell. The CpG island located at the promoter region of BP1 (CpG island 1) was almost unmethylated in a representative CML patient and K562 cell as detected by bisulfite sequencing. A white circle indicates unmethylated CpG dinucleotide, whereas a black circle indicates methylated CpG dinucleotide. Each line represents an independent clone-sequencing result of BSP product of K562 cell or CML patient. d The discriminating value of BP1 in CML patients. BP1 expression may serve as a potential biomarker for distinguishing CML patients from controls with an AUC value of 0.624. e The discriminating value of DLX7 in CML patients. DLX7 expression may serve as a potential biomarker for distinguishing CML patients from controls with an AUC value of 0.699 stages (P = 0.008,

DLX4 isoforms BP1 and DLX7 showed opposite function in cell division
Since the DLX4 isoforms BP1 and DLX7 showed a differentially expressed pattern in CML, we next investigated the potential role of these two DLX4 isoforms in cell division. Because it is hard to knockdown BP1 and DLX7 expression alone by RNA interference, we performed gain-of-function experiments in K562 cells with BP1/DLX7 overexpression by infection with lentivirus (Genechem, Shanghai, China) carrying the complete CDS of BP1 and DLX7. The establishment of K562 cells stably overexpressing BP1 was proved by both RT-qPCR (Fig. 2a, b) and western blot (Additional file 12), whereas K562 cells stably overexpressing DLX7 were only confirmed by RT-qPCR because of the lack of specific antibodies (Fig. 2c, d). As expected, overexpression of BP1 significantly increased the proliferation rate of K562 cells together with S/G2 promotion (Fig. 2f-h), whereas DLX7 overexpression markedly decreased the proliferation rate of K562 cells together with G1 arrest (Fig. 2i-k). Collectively, the functional experiments in vitro demonstrated the promitotic effects of BP1 but the antimitotic effects of DLX7 in leukemogenesis.
(See figure on next page.) Fig. 2 DLX4 isoforms BP1 and DLX7 showed opposite functions in leukemogenesis both in vivo and in vitro. a BP1 expression in K562 cells before and after BP1 transfection. BP1 expression was significantly increased after BP1 transfection in K562 cells as detected by real-time quantitative PCR (RT-qPCR). K562-NC (control) was defined as 100% expression for BP1 transcript. b DLX7 expression in K562 cells before and after BP1 transfection. DLX7 expression was not changed after BP1 transfection in K562 cells as detected by RT-qPCR. K562-NC (control) was defined as 100% expression for DLX7 transcript. c DLX7 expression in K562 cells before and after DLX7 transfection. DLX7 expression was significantly increased after DLX7 transfection in K562 cells as detected by RT-qPCR. K562-NC (control) was defined as 100% expression for DLX7 transcript. d BP1 expression in K562 cells before and after DLX7 transfection. BP1 expression was not changed after DLX7 transfection in K562 cells as detected by RT-qPCR. K562-NC (control) was defined as 100% expression for BP1 transcript. e The proliferation ability in K562 affected by BP1 overexpression. K562 cells with BP1 overexpression (K562-BP1) showed a significantly increased proliferation rate than those without BP1 overexpression (K562-NC). f The cell cycle in K562 affected by BP1 overexpression. K562-BP1 cells showed a significantly higher percentage of S/G2 phase than K562-NC cells. g Representative cell cycle diagram of K562-NC and K562-BP1 cells by flow cytometry. h The proliferation ability in K562 affected by DLX7 overexpression. K562 cells with DLX7 overexpression (K562-DLX7) showed a markedly reduced proliferation rate than those without DLX7 overexpression (K562-NC). i Cell cycle in K562 affected by DLX7 overexpression. K562-DLX7 cells showed significantly lower percentage of S/G2 phase than K562-NC cells. j Representative cell cycle diagram of K562-NC and K562-DLX7 cells by flow cytometry. k The flow chart of the in vivo experiment. l The tumor load in NCG mice affected by K562 cells with BP1 and DLX7 overexpression. The tumor load of K562-BP1 group mice was significantly higher, whereas the tumor load of K562-DLX7 group mice was significantly lower compared with the tumor load of K562-NC group mice via bioluminescence imaging at the 28th day and 42th day. m The representative tumor load diagram of NCG mice with K562-NC, K562-BP1, and K562-DLX7 cells injection as detected by bioluminescence imaging. n The representative tumor volume of NCG mice with injection of K562-NC, K562-BP1, and K562-DLX7 cells. *P < 0.05; **P < 0.01; ***P < 0.001; NS, no significance

Functional role of DLX4 isoforms BP1 and DLX7 in a mouse xenograft model
To further determine the role of DLX4 isoforms BP1 and DLX7 in vivo, we constructed xenograft mouse models by tail vein injection of K562-NC/K562-BP1/K562-DLX7 cells in NCG mice. The experimental procedures are shown in Fig. 2m. Firstly, the progression of the tumor load in mice was tracked using bioluminescence imaging, revealing that the tumor load of the K562-BP1 group mice was higher, whereas the tumor load of K562-DLX7 group mice was lower as compared with the tumor load of K562-NC group mice (Fig. 2n, o). Moreover, the tumor size in K562-BP1/ K562-DLX7/K562-NC group mice is shown in Fig. 2p. Overall, the functional studies in vivo further support the tumor-promoting role of BP1 but the tumor-suppressing role of DLX7 in leukemogenesis.

Molecular mechanism of DLX4 isoforms BP1 and DLX7 in leukemogenesis
Firstly, to explore how the transcription factor BP1 exerts its function, we performed ChIP-Seq and RNA-Seq to identify candidate downstream genes in leukemogenesis (Fig. 3a). Using RNA-Seq, a total of 2834 genes were identified to be differentially expressed after BP1 overexpression (Fig. 3b and Additional file 2: Table S2, Additional file 3: Table S3 and Additional file 4: Table S4). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the differentially expressed genes (DEGs) were significantly enriched in several cancer-related signaling pathways such as MAPK and AMPK signaling (Fig. 3c). Moreover, ChIP-Seq results revealed that the transcription factor BP1 may bind to the promoter region of 661 genes (Additional file 5: Table S5). Motif analysis identified the top three motifs in BP1-bound promoters, which is shown visually in Fig. 3d. Venn analysis of the RNA-Seq and ChIP-Seq results showed that a total of 91 genes including 79 mRNA and 12 lncRNAs may act as direct downstream targets of BP1 (Fig. 3e, f ). Four mRNAs (RREB1, VEGFA, ASAP1, and NPHP4) as well as three lncRNAs (ID2-AS1 transcript 1, ID2-AS1 transcript 2, and SGMS1-AS1) were further validated by RT-qPCR (Fig. 3g). To further confirm the biological association of RREB1 and SGMS1-AS1 with BP1, we next performed rescue experiments during luekemogenesis caused by BP1 overexpression. As expected, knockdown of RREB1 and  -AS1 expression by siRNA markedly impaired the proliferation after BP1 overexpression in K562 cells (Fig. 3h).
Secondly, because of a lack of specific antibodies, we used ATAC-Seq and RNA-Seq to identify candidate downstream genes in leukemogenesis (Fig. 3i). RNA-Seq revealed a total of 2277 genes identified as differentially expressed after DLX7 overexpression ( Fig. 3j and Additional file 6: Table S6, Additional file 7: Table S7, Additional file 8: Table S8). KEGG analysis revealed that the DEGs were significantly enriched in several cancer-related signaling pathways such as FoxO signaling (Fig. 3k). Additionally, DLX7 may bind to the promoter region of 9361 genes by ATAC-Seq (Additional file 9: Table S9). Motif analysis identified the top one motif in DLX7-bound promoters, as shown visually in Fig. 3l. Venn analysis of the RNA-Seq and ATAC-Seq results showed that a total of 282 genes including 151 mRNA and 131 lncRNAs may act as direct downstream targets of DLX7 (Fig. 3m, n). Four mRNAs (PTPRB, PDZK1, DLC1, and CASP9) as well as two lncRNAs (NEAT1 and LINC-PINT) were further validated by RT-qCR (Fig. 3o). To further confirm the biological connections of PTPRB and NEAT1 with DLX7, we next performed rescue experiments during luekemogenesis caused by DLX7 overexpression. As expected, knockdown of PTPRB and NEAT1 expression by siRNA partially revived the proliferation after DLX7 overexpression in K562 cells (Fig. 3p).

Identification of an lncRNA-miRNA-mRNA network in leukemogenesis caused by BP1
overexpression Given the results above, most of the differentially expressed mRNAs together with miRNAs that were not identified as direct downstream targets of BP1 greatly aroused our attention. Based on the fact that several lncRNAs were verified as downstream targets of BP1, we deduced that the miRNAs and mRNAs were downstream targets of the activated lncRNAs and may be indirectly activated by BP1. It is well known that lncR-NAs located in cytoplasm may work through a competing endogenous RNA (ceRNA) (lncRNA-miRNA-mRNA) network. Herein, we verified that an lncRNA SGMS1-AS1 directly activated by BP1 was anticipated to be located mainly in the cytoplasm of all available cell types by using lncLocator (http:// www. csbio. sjtu. edu. cn/ bioinf/ lncLo cator/) (Fig. 4a). Further LncRNA-FISH assays confirmed that SGMS1-AS1 was localized mainly in the cytoplasm of K562 cells (Fig. 4b). To explore the downstream miR-NAs, we identified miRNAs that can bind with SGMS1-AS1 by using publicly available online tools (LncBase Predicted v.2 and starBase) together with miRNA sequencing. The results revealed that only miR-181d-5p was shared in all three conditions (Fig. 4c and Additional file 10: Table S10). Additionally, the reduced expression of miR-181d-5p after BP1 overexpression in K562 cells was further validated by RT-qPCR (Fig. 4d). Moreover, the dual-luciferase reporter assays revealed that overexpression of miR-181d-5p significantly reduced the luciferase activity of the wild-type SGMS1-AS1 vector but not of the mutated SGMS1-AS1 vector (Fig. 4e). Together, these findings prove a direct association between SGMS1-AS1 and miR-181d-5p in leukemogenesis.
Because Ago2 is a core component of the RNA-induced silencing complex (RISC) involved in miRNA-mediating mRNA destabilization or translational repression, we performed RIP assays by using an anti-Ago2 antibody, which demonstrated that endogenous SGMS1-AS1 together with miR-181d-5p and SRPK2 was preferentially enriched in Ago2-RIPs compared with control IgG-RIPs (Fig. 4i, j). All these findings suggest that SGMS1-AS1 functions via ceRNA network to play its biologic role during leukemogenesis.
Combining these findings, a schematic diagram of the roles of DLX4 gene isoforms in promoting leukemogenesis is presented in Fig. 5.

Discussion
The expression pattern of DLX4 and its clinical significance have been demonstrated in several human cancers. Gao et al. revealed that DLX4 expression was dramatically increased in patients with hepatocellular carcinoma and was associated with poor prognosis [16]. Man et al. demonstrated that overexpression of BP1 independently affected disease-free survival in patients with NSCLC [17]. Moreover, mounting evidence has shown the adverse effect of DLX4/BP1 overexpression on clinical outcome in patients with breast cancer [12][13][14][15]. Importantly, Haga et al. also reported the overexpression pattern of BP1 in all types of acute leukemia [18]. However, several studies have also shown the hypermethylation of DLX4 with a role of transcriptional inactivation in stage I NSCLC, uterine cervical LSILs, breast cancer, and CLL [19][20][21][22][23]. Moreover, DLX4 hypermethylation was found to be associated with disease progression in uterine cervical LSILs and NSCLC [22,23]. Taking these results together, the expression pattern of DLX4 seems to be "contradictory" when contrasted with the methylation condition of DLX4 in human cancers, especially breast cancer. In-depth analysis in previous studies regarding DLX4 methylation has focused on the CpG islands located in the promoter region of DLX7 but not the promoter region of BP1 [19][20][21][22][23]. The methylation pattern of the CpG island located at the promoter region of BP1 has rarely been studied. Our preliminary studies also confirmed the phenomenon of DLX4 hypermethylation in MDS, AML, and CML [24][25][26]. Moreover, DLX4 methylation with its role in silencing DLX7 expression but not BP1 expression was further verified in AML and CML [25,26]. In addition, BP1 and DLX7 showed opposite expression patterns in AML patients; that is, BP1 expression was increased, whereas DLX7 expression was decreased [27]. In the work described herein, we further detected the expression pattern of BP1 and DLX7 in CML patients, with results similar to AML patients [27]. Moreover, the CpG islands located at the promoter region of BP1 were nearly unmethylated in CML. These results further support that the DNA methylation-mediated differential expression pattern of DLX4 isoforms may play a different role in leukemogenesis.
Regarding the role of DLX4 in tumorigenesis, most previous studies have mainly focused on DLX4 isoform 1, and thus BP1 is also called DLX4 in some studies. Zhang et al. demonstrated that DLX4 activated expression of TWIST to promote epithelial-tomesenchymal transition (EMT), cancer migration, invasion, and metastasis in breast cancer [33]. Alternatively, BP1 overexpression markedly enhanced cell proliferation and metastatic potential in estrogen receptor (ER)-negative Hs578T breast cancer cells [34]. Haria et al. showed that DLX4 induced CD44 expression by stimulating interleukin (IL)-1β-mediated nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) activity, thereby promoting peritoneal metastasis of ovarian cancer [35]. Moreover, DLX4 played an oncogenic role in clear cell renal cell carcinoma via inducing proliferation and EMT, and was associated with poor prognosis [36]. DLX4 isoform BP1 overexpression promoted cell proliferation, migration in endometrial cancer with clinically prognostic effect [37]. However, the direct role of DLX7 was poorly determined. In the current work, we systemically investigated the direct role of BP1 and DLX7 during leukemogenesis by in vivo and in vitro studies. Our study demonstrated that DLX4 isoforms BP1 and DLX7 have opposite functions in leukemogenesis, with BP1 playing an oncogenic role whereas DLX7 plays a tumor suppressive role in leukemogenesis. Interestingly, DLX4 isoforms BP1 and DLX7 have distinct functions in the regulation of the b-globin gene [38].
The potential mechanism of DLX4, as a transcription factor, in carcinogenesis has been investigated. Most previous studies have explored the potential downstream targets of DLX4 by single-gene identification based on ChIP-PCR. For instance, Kluk et al. disclosed that BP1 homeoprotein repressed BRCA1 expression by direct binding to its first intron in sporadic breast cancer [39]. Alternatively, BP1 transcriptionally activates bcl-2 and inhibits tumor necrosis factor (TNF)-α-induced cell death in MCF7 breast cancer cells [40]. Trinh et al. showed that DLX4 blocked the antiproliferative effect of transforming growth factor (TGF)-β by inhibiting TGF-β-mediated induction of p15 Ink4B and p21 WAF1/Cip1 expression and inducing expression of c-myc independently of TGF-β/Smad signaling [41]. Moreover, DLX4 promoted nasopharyngeal carcinoma progression via inducing YB-1 expression [42]. In this study, we used ChIP-Seq/ATAC-Seq together with RNA-Seq to identify the underlying targets of DLX4 in leukemogenesis by whole-genome-scale screening. Interestingly, although the two isoforms BP1 and DLX7 are predicted to share a common homeodomain, the binding motifs were discrepant for BP1 and DLX7. A possible mechanism is that the locus occupancy is mediated by other proteins that bind DNA directly. Berger et al. revealed that variation in homeodomain DNA binding sequence recognition may be a factor in functional diversity and evolutionary success [43]. Similarly, Song et al. using ChIP on chip (ChIP-on-chip) and gene expression microarray assays to screen the downstream targets of BP1 in ER-breast cancer cells [44]. A total of 18 genes were identified and verified, of which some were involved in a variety of tumorigenic pathways in breast cancer development and progression [44]. Although VEGFA was identified as a downstream target of DLX4 in the studies by Song et al. [44] and Hara et al. [45] as well as our study, most of the genes were not identified by our investigation. These results suggest that BP1 could regulate different downstream target among different types of human cancer.
On the basis of whole-genome study, most of the differentially expressed mRNAs together with miRNAs that were not identified as direct downstream targets of BP1 greatly aroused our attention. An emerging phenomenon is that one of the molecular mechanisms by which lncRNAs regulate gene expression is to interact with miRNA as ceRNAs that bind to miRNA response elements and protect miRNAs from binding to and repressing target mRNAs [46][47][48]. Accordingly, we deduced that the miRNAs and mRNAs were downstream targets of the activated lncRNAs and may be indirectly activated by BP1. In this study, we identified a BP1/SGMS1-AS1/miR-181d-5p/SRPK2 ceRNA network in leukemogenesis. Although there are no studies regarding SGMS1-AS1 in leukemogenesis, several studies have indicated the role of miR-181d and SRPK2 in AML biology. Su et al. reported the direct role of the miR-181 family in normal hematopoiesis and AML development [49]. Moreover, miR-181d/RBP2/NF-κB p65 feedback regulation promoted CML progression into blast crisis [50]. Additionally, Jang et al. revealed that SRPK2 promoted leukemia cell proliferation by phosphorylating acinus and regulating cyclin A1 [51]. Collectively, all these results suggest the crucial role of the BP1/SGMS1-AS1/miR-181d-5p/SRPK2 network in leukemogenesis.

Conclusion
Taken together, our findings reveal that DNA methylation-mediated differential expression of DLX4 isoforms BP1 and DLX7 have opposite functions in leukemogenesis. BP1 plays an oncogenic role in leukemia development, whereas DLX7 acts as a tumor suppressor gene. Moreover, the SGMS1-AS1/miR-181d-5p/SRPK2 ceRNA network activated by BP1 plays a significant role in leukemogenesis. These results suggest DLX4 as a therapeutic target in antileukemia therapy.