We recruited 23 MDS patients at the Department of Hematology, Huashan Hospital, Fudan University (Shanghai, China) from January 1, 2015 to December 30, 2016. The inclusion criteria consisted of newly diagnosed patients and no previous treatment and excluded patients with myelodysplastic/myeloproliferative neoplasms (MDS/MPN) and MDS which had progressed to acute myeloid leukemia. Bone marrow samples were collected for all patients and patients were diagnosed with MDS and classified according to the 2008 revision of the World Health Organization (WHO) criteria (22). There were 15 men and 8 women with a median age of 65 years (range, 29–82 years), in the MDS group. In addition, bone marrow samples from 7 cases with non-hematological malignancies with a median age of 38 years (range 23–78 years) were obtained between January 1, 2015 and December 30, 2015, and were analyzed as controls. No statistical difference was found with regards to age and sex distributions among the subjects with MDS and non-MDS controls. For microarray analysis groups, we obtained bone marrow samples from anther 4 MDS patients (2 patients with refractory anemia with excess blasts 1 (RAEB-1) and 2 patients with RAEB-2] and 4 age-matched patients with hypersplenism from the Department of Hematology, Huashan Hospital, Fudan University (Shanghai, China) in 2014. We collected 2 ml of all bone marrow samples using heparin and stored the samples at room temperature within 8 h. Total RNA was isolated from bone marrow samples within 8 h and then stored at −80°C. The present study was approved by the Ethics Committee of Huashan Hospital, School of Medicine, Fudan University. All participants provided a written informed consent form before being enrolled into this study.

Total RNA was isolated from bone marrow samples using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.). The purity of RNA was assessed by measuring the optical density (OD) 260/280 ratio using a NanoDrop spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.). The RNA samples were reverse transcribed into cDNA using Takara PrimeScript RT Master Mix (Takara Bio Inc., Otsu, Shiga, Japan). The 10-µl reaction mixture comprised 2 µl 5X PrimeScript RT Master Mix, 500 ng RNA and RNase-Free distilled H₂O, which was used to ensure that the total reached 10 µl. Subsequently, the mixture was amplified at 37°C for 15 min, 85°C for 5 sec and was stored at 4°C for further use.

For RT-qPCR, cDNA samples were amplified using an Applied Biosystems™ 7500 Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.) in a 20-µl mixture containing 10 µl SYBR^® Premix Ex Taq™ (Takara Bio, Inc., Otsu, Japan), 0.4 µl of each primer (10 µmol/l), 0.4 µl ROX Reference Dye II (Takara Bio, Inc.), 2 µl cDNA, and 6.8 µl ddH₂O at 95°C for 30 sec, followed by 40 cycles at 95°C for 5 sec and 60°C for 34 sec. Gene expression levels were expressed relative to the expression of GADPH. The primers were designed with the Primer 3.0 online software (http://bioinfo.ut.ee/primer3-0.4.0/) and synthesized by BioTNT (Shanghai, China). The primer sequences were as follows: lncENST00000444102, 5′-TATCGACGTAGTTAAAGCCCACT-3′ and 5′-CTTCTGCCCTTCACATCCTCT-3′; ABAT, 5′-CCGACTACAGCATCCTCTCC-3′ and 5′-GGTTCTCTTTCACAAACTCTTCC-3′; GAPDH, 5′-GACCTGACCTGCCGTCTA-3′ and 5′-AGGAGTGGGTGTCGCTGT-3′. The relative levels of gene expression were calculated using the 2^−ΔΔCq method (23).

Subsequently, cDNA samples from microarray analysis groups were used to profile differentially expressed gene transcripts and lncRNAs using Agilent human genome-wide gene expression BeadChips (Agilent Technologies, Inc., Santa Clara, CA, USA) according to the manufacturer's protocol at 65°C. The expression values were normalized with Robust Multi-array Average (RMA) of background-adjusted, normalized and log₂-transformed using the statistical software package R (24).

After profiling each DEL and DEM in the MDS samples, we specifically identified the altered ABAT expression in MDS patients for construction of the ABAT-DEL-DEM network. In brief, the co-expression of ABAT and particular DELs and lncRNA-mRNA correlation were evaluated using the Pearson correlation coefficient (PCC). The ABAT-DEL-DEM co-expression network was further filtered by the overlapping DELs using software of Cytoscape 3.4.0 (The Cytoscape Consortium, San Diego, CA, USA. Weblink: http://cytoscape.org/download.html) according to a previous study (25).

Subsequently, we investigated the potential role of the identified ABAT-DEL-DEM co-expression network using the GO and KEGG pathway and network analyses using the Database for Annotation, Visualization and Integrated Discovery (DAVID, v6.8 tool) (26,27). The KEGG pathway enrichment analysis was performed using the clusterProfiler package in R/bioconductor. Furthermore, the key KEGG signaling pathway was analyzed using the R package pathview (28).

For the bioinformatic analysis of lncRNAs, we first accessed the UCSC database (http://genome.ucsc.edu) to prioritize ABAT-associated lncRNAs and searched the NCBI human genomes database to identify the chromosomal localization of each lncRNA. We also utilized PhyloCSF, a novel comparative genomics method to analyze multispecies nucleotide sequence alignment to determine whether DNA sequences where each lncRNA resides are likely to represent a conserved protein-coding or non-coding region (29). PhyloCSF scores for selected phylogenies may be viewed in the UCSC Genome Browser by copying the URL ‘http://www.broadinstitute.org/compbio1/PhyloCSFtracks/trackHub/hub.txt’ into the ‘My Hubs’ tab under ‘track hubs’. PhyloCSF outputs a score, positive if the alignment is likely to represent a conserved coding region, and negative otherwise. Moreover, we also used computational and mathematical methods to predict advanced lncRNA structure (http://rna.tbi.univie.ac.at/RNAfold/UyEBF8akiV), and used the TRANSFAC database to predict the transcription factor binding sites (TFBS) of each lncRNA (http://www.gene-regulation.com) (30).

A human AML cell line (THP-1) was purchased from Chinese Academy of Sciences (Shanghai, China) and a human AML cell line transformed from MDS cells (SKM-1) was obtained from the Japanese Collection of Research Bioresources (LCRB; Tokyo, Japan). 293T cells were used as a control (Chinese Academy of Sciences, Shanghai, China). All cell lines were cultured in Roswell Park Memorial Institute-1640 medium (RPMI-1640; Hyclone; GE Healthcare Life Sciences) supplemented with 10% fetal bovine serum (Gibco-BRL; Thermo Fisher Scientific, Inc.) in a humidified incubator with 5% CO₂ at 37°C.

To confirm our ABAT-DEL-DEM co-expression network, we constructed lentiviral vectors and prepared lentivirus to stably overexpress lncRNA or knock down ABAT expression. All plasmid vectors were produced using standard cloning techniques (31). The lncRNA was overexpressed and cloned into GV470 lentiviral vectors (GeneChen, Shanghai, China). The shRNA hairpins that targeted the 3′-untranslated region (3′-UTR) of ABAT were designed and cloned into pGMLV-SC5 lentivirus vectors (Genomeditech, Shanghai, China). The sequences of oligonucleotides were 5′-GCTGGAGACGTGCATGATTAA-3′. The SKM-1 and THP-1 cells were respectively cultured in 6-well plates at a density of 1×10⁶ cells/ml overnight and transduced with lncRNA overexpression lentivirus and shRNA lentivirus plus a scramble lentivirus (negative controls), and cells without any lentivirus transduction were considered as controls. Since the lentivirus had green fluorescence, we evaluated transduction efficiency by flow cytometry at 72 h. The percentage of positive cells was >80% by flow cytometry. The growth medium was then replaced with 1 or 2 µg/ml puromycin, respectively, to select stable cells for two weeks. These cells were then subjected to quantitative reverse transcription polymerase chain reaction (RT-qPCR) to verify expression levels of lncRNA and ABAT gene.

Cell viability was assessed using the Cell Counting Kit-8 (CCK-8) assay kit from Dojindo Laboratories (Kumamoto, Japan). In brief, stably transfected and control cells were seeded into a 96-well culture plate at a density of 1×10³ cells/well and cultured for up to five days, with media replaced every three days. At the end of each experiment, 5 µl of the CCK-8 reagent was added into each well and the cells were further cultured for 4 h. Then, the optical density of cells was measured at 450 nm using a spectrophotometer (Tecan, Männedorf, Switzerland). The experiments were performed in triplicate and repeated at least three times.

Cells cultured for 48 h (6-well plate, 1×10⁶ cells/ml) were collected and washed twice with phosphate-buffered saline (PBS), and then re-suspended in 200 µl of the binding buffer. Annexin V-APC (5 µl) and 7-AAD (5 µl; Thermo Fisher Scientific, Inc.) were added and the mixture was incubated in the dark at room temperature for 15 min. The rate of cellular apoptosis was then measured using a flow cytometer (BD Accuri C6; BD Biosciences, San Jose, CA, USA). Data were statistically analyzed using the software Flowjo 7.6 (Flow Jo, LLC, Ashland, OR, USA). The experiments were performed in duplicate and repeated at least once.

Student's t-test was used to identify DELs and DEMs between patients with MDS and hypersplenism by calculation of P<0.05 and fold change (FC)>2 for each DEM and DEL. The correlation of ABAT-DEL-DEM co-expression was evaluated using Pearson correlation coefficient (PCC). Pairs with a PCC threshold >0.95 and P-value <0.05 were used to construct the ABAT and DEL co-expression association matrix. Pairs with a PCC threshold >0.99 and P-value <0.01 were selected as the meaningful value to construct the DEL and DEM co-expression network. A P-value <0.05 using Fisher's exact test and a kappa (κ) coefficient of 0.4 were used as threshold values in GO and KEGG pathway and network analyses. The in vitro experimental data were analyzed using Graphpad Prism 6 software (GraphPad Software, La Jolla, CA, USA). Comparisons between two groups were analyzed by Student's t-test when data conformed to normal distribution, if not, a non-parametric Kruskal-Wallis test was performed. Following a Kruskal-Wallis test, Dunn's Multiple Comparisons test was used as the post test to compare the difference in the sum of ranks between two columns with the expected average difference. A P<0.05 was considered statistically significant.

The microarray data were then scanned using SureScan Scanner (Agilent Technologies) and analyzed as previously reported (8). Bone marrow cells from 4 MDS patients were compared with cells from 4 age-matched hypersplenism controls, and the microarray analysis revealed a total of 543 DELs and 2,705 DEMs. Among them, 285 (52.5%) DELs were downregulated and 258 (47.5%) DELs were upregulated in MDS patients, whereas 1,521 (56.2%) DEMs were downregulated and 1,184 (43.70%) DEMs were upregulated in MDS. The volcano plots and heatmaps provide an overview of the DEL and DEM microarray data, respectively (Fig. 1).

The ABAT-DEL-DEM co-expression network was constructed using the Cytoscape program. To generate the ABAT-DEL-DEM co-expression network, we related ABAT expression in MDS samples to each of the DELs using the PCC test. We were able to select six DELs (Table I) that were co-expressed with ABAT with a PCC threshold >0.95 and P<0.05 and 135 co-expressed mRNAs with a PCC threshold >0.99 and P<0.01 in MDS to construct the lncRNA and mRNA co-expression network (Table SI). Next, their potential interaction was determined using the Cytoscape program. The data showed that the co-expression network was composed of six DELs related to ABAT and 135 co-expressed mRNAs (Fig. 2). In this network, three lncRNAs were up-regulated (lncRNA1, lncRNA2 and lncRNA6), whereas the other three were downregulated (lncRNA3, lncRNA4 and lncRNA5). The network showed that a particular mRNA could correlate with numerous lncRNAs, while a single lncRNA was also able to correlate to various mRNAs, implying that an inter-regulation of lncRNAs and mRNAs occurs in MDS.

We next performed analyses using GO terms and the KEGG pathways and networks, and found that the co-expression network was mostly enriched in several biological processes, such as cellular components and molecular functions. In particular, the GO analysis revealed that the co-expression network mainly participated in response to organic cyclic compounds, cell proliferation, cell part morphogenesis, regulation of cell proliferation, and enzyme-linked receptor protein signaling pathway (Fig. 3A) for molecular functions of protein homodimerization activity. The GO term analysis also revealed that the co-expression network is localized in different areas within cells, such as the lysosome, vacuole, cell projection, and extracellular vesicular exosome. The KEGG data showed that the co-expression network was involved in different pathways, such as the phagosome and metabolic pathways (Fig. 3B). The differentially expressed gene of neutrophil cytosolic factor 2 (NCF2) plays role in neutrophil phagosome. In addition, the genes of serine palmitoyltransferase long chain base subunit 2 (SPTLC2), 4-aminobutyrate aminotransferase (ABAT) and quinolinate phosphoribosyltransferase (QPRT) were key metabolic enzymes which play a role in metabolic pathways (https://www.ncbi.nlm.nih.gov/gene) (Table II).

In this study, we revealed six lncRNAs and ranked lncENST00000444102 as the most important lncRNA with the highest correlation index related to the ABAT gene (Table I). lncENST00000444102 (lncRNA4) was co-expressed with the genes, NCF2, SPTLC2, ABAT and QPRT which all take part in the pathways (Table II and Fig. 2), thus we focused on lncENST00000444102 for further research. We then used tools in the UCSC database (http://genome.ucsc.edu) to re-annotate them and focused on chr6: 167382710-167411729, reverse strand (Fig. 4A). This lncRNA (ENST00000444102, named lncENST00000444102) was one of the most downregulated lncRNAs in MDS, residing on chromosome 6 in the human genome as an antisense RNA, consisting of three exons and spanning nearly 29 kilobases (kb). The advanced structure of lncENST00000444102 was predicted by RNAfold WebServer. We typed in the RNA sequence and chose the options of fold algorithms, and the minimum free energy (MFE) structure and the centroid secondary structure were obtained (Fig. 4B) (32). PhyloCSF scores for lncENST00000444102 were viewed in the UCSC Genome Browser tracks. We found that the scores were negative, with which we confirmed that the lncENST00000444102 was a non-coding RNA (Fig. 5). The gene was located among chemokine receptor 6 (CCR6), FGFR1 oncogene partner (FGFR1OP), and ribonuclease T2 (RNASET2) (Fig. 5), indicating that this particular lncRNA may be able to regulate the transcription of these three genes in cells. To confirm this hypothesis, we predicted the transcription factor binding sites (TFBS) of lncENST00000444102 using tools listed in the TRANSFAC database and found that lncENST00000444102 was able to target chicken ovalbumin upstream promoter-transcription factor (COUP-TF), which plays an important role in the regulation of organogenesis, neurogenesis, metabolic homeostasis, and cellular differentiation during embryonic development, and hepatocyte nuclear factor 4 (HNF-4), mainly playing a role in hepatic diseases (Table III).

The alterations of these two molecules were next confirmed in MDS bone marrow samples from 23 MDS patients, 7 non-hematological malignancies, and two cell lines (Fig. 6). lncENST00000444102 expression was found to be significantly downregulated in MDS patients (P<0.0001, Fig. 6A) when compared to non-MDS cases as well as in the SKM-1 and THP-1 cell lines (P<0.0001, Fig. 6B) when compared to HEK-293 cells. ABAT expression was also downregulated in the bone marrow from MDS cases (P<0.001, Fig. 6C). We then overexpressed lncENST00000444102 and knocked down ABAT expression in SKM-1 and THP-1 cells (Fig. 7A and B). We found that stable overexpression of lncENST00000444102 induced ABAT expression by 2-fold in the SKM-1 and THP-1 cell lines (P<0.01, Fig. 7C). However, when ABAT expression was stably knocked down in the SKM-1 and THP-1 cell lines, the level of lncENST00000444102 was not significantly (ns) altered (P>0.05, Fig. 7D).

The cell viability CCK-8 assay showed that after 24 h in culture, the viability of SKM-1 and THP-1 cells with stable lncENST00000444102 overexpression started to decrease when compared with that of the control (P<0.05, Fig. 8A). The fraction of apoptotic cells was 22.41±2.596 in the lncENST00000444102-overexpressing SKM-1 cells, and 8.650±0.889 in the negative control; the fraction of apoptotic cells was 20.58±2.190 in the lncENST00000444102-overexpressing THP-1 cells and 8.192±0.997 in the negative control group (P<0.001, Fig. 8B). The flow cytometric apoptosis assay revealed that lncENST00000444102 overexpression promoted tumor cells to undergo apoptosis compared to control cells (P<0.001, Fig. 9). These data revealed that lncENST00000444102 overexpression reduced MDS cell viability by induction of MDS cell apoptosis.