Establishing an ir-lncRNAs pairs for AML patients
The expression profiles of 1289 ir-lncRNAs were used for differential expression analysis (Fig. 1A), and there were total 380 differentially expressed genes in AML and normal tissues, of which 210 genes were upregulated and 170 were downregulated (Fig. 1B). Subsequently, we constructed a 0 or 1 matrix to generate ir-lncRNAs pairs based on above 380 genes and 16381 pairs were identified using univariate analysis, and 1710 ir-lncRNAs pairs were verified by Lasso penalized Cox regression to screen the prognosis-related genes with potentiality. As long as genes with non-zero coefficients had prognostic value in the Lasso penalized regression model, and eventually, we selected one model which produced a group of eight lncRNAs pairs (AC093642.1|HOXA10.AS, PCED1B.AS1|AL365361.1, AL158210.1|PCAT18, U62631.1|HOXA10.AS, AL365361.1|LINC00205, AL137779.1|AC079305.1, ILF3.DT|MIR222HG, AC018690.1|AC022210.1)(Figs. 1C, D). To validate the Lasso penalized Cox regression model, univariate and multivariate Cox proportional hazard regression analyses were performed (Fig. 1E, F).
Construction of the ir-lncRNAs pairs prognostic model. (A) Heatmap of the expression of lncRNAs. (B) Differentially expressed genes in AML and normal tissues. (C) The LASSO coefficient values at various levels of penalty. (D)The confirmation of the best lambda value by LASSO Cox regression analysis. (E, F) Univariate and multivariate Cox regression analysis of the each pairs in this model.
Evaluation Of The Eight Ir-lncrnas Pairs Risk Model
The AUCs of 1-year ROC curve for the ir-lncRNAs pairs was 0.825 and the optimal cutoff value was 2.201 (Fig. 2A). 186 patients in the training dataset were divided into high-risk and low-risk groups based on the cutoff point, and the K-M survival analysis presented a much worse outcome in the high-risk group than that of the low-risk group (log-rank test, p < 0.0001, Fig. 2B). The risk score and survival data of each patient were shown in Fig. 2C, the results confirmed that death was more frequently observed in the high-risk group set than the low-risk group. And the AUCs of a time-dependent ROC curve of 1-, 2- and 3-year calculated by the risk score model were 0.825, 0.883 and 0.928 respectively (Fig. 2D).
The testing column (TCGA, 130 patients) was used to verified the predictive values of the eight ir-lncRNAs pairs risk score. Patients were also divided into low-risk and high-risk groups based on the cutoff value and K–M curve promoted patients classified to these two level risk score obtained an obviously different OS (log-rank test, p = 0.0078, Fig. 2E), and the risk score and survival data of each patient are shown in Fig. 2F. The AUCs of 1-, 2- and 3-year ROC curves were 0.657, 0.632 and 0.677 (Fig. 2G), these results showed that the ir-lncRNAs pairs prognostic model might be a potential predictor to judge the OS of patients with AML.
Evaluation of the ir-lncRNAs pairs prognostic model in training and testing column. (A) AUC of 1-year ROC curve and confirmation of cutoff value. (B) Time-dependent ROC curve analysis for the prediction survival using the risk model in training dataset. (C) Dot plots comparing the outcomes of subjects and survival time in the high- and low-risk groups. (D) Survival curves of two groups were plotted using the K-M method. (E-G) Results of testing dataset.
Relationship Between Clinical Factors And Ir-lncrnas Pairs Risk Model
The strip chart showed that the distribution of CEBPA mutation (p < 0.05) and ELN risk stratification (p < 0.001) were distinct in high- and low-risk patients (Fig. 3A). Meanwhile, we divided AML patients into several subtypes according to the gender, European leukemia network (ELN) risk stratification, mutational status of the FLT3-ITD, NPM1 and CEBPA, respectively, then compared the risk score of different groups. The ladder diagrams showed that CEBPA mutated group had lower risk score (Fig. 3E), and similarly as expected, the difference in risk score between patients with favorable, intermediate and adverse groups were statistically significant (Fig. 3F). There were no difference across different gender, NPM1 or FLT3-ITD mutations(Figs. 3B, C, D). In subgroup analyses, no clear separation in OS K-M curves was observed between male and female subgroups (Figs. 3G, H) or between NMP1 and FLT3 wild-type and mutated AML subgroups (Figs. 3I-L).
Relationship between clinical factors and ir-lncRNAs pairs risk model. (A) Strip chart showing that distribution of CEBPA mutation and ELN risk stratification are significantly different in high- and low-risk groups. (B-F) Scatter diagrams comparing the risk score of each patients in different clinical features. (G, H) K–M survival curves showing the subgroup analyses of gender, NMP1 (I, J) and FLT3 (K, L) mutation status.
Differences In Tumor Immune Landscape Between High- And Low-risk Patients
We analyzed the distribution of mutations across the different genes differed between high- and low-risk groups in the Vizome database. Mutation frequencies per molecular in different groups were show in Figs. 4A-C. And according to statistical results, there were 384 genes mutation occurred only in low-risk groups, 668 genes only in high-risk patients, and 71 genes in both groups (Fig. 4D).
In the meantime, we assessed the difference in expression of 34 HLA-related genes in all three datasets between high- and low-risk patients and found that HLA-B, HLA-E, HLA-F, HLA-T, HLA-U, HLA-DOB, HLA-DQA1 and HLA-DQB1 were elevated (P < 0.05) in the high-risk signature group (Fig. 4E). In addition, the expression of 37 immune checkpoint molecules were performed between high- and low-risk groups and we found 14 molecules had a higher expression in high-risk patients, including CTLA-4, PDCD1, CD274, LAG3 and PDCD1LG2, etc (Fig. 4F), indicating that high-risk patients may respond better and had better outcome when receiving immune checkpoint inhibitors.
Differences in tumor immune landscape between high- and low-Risk Patients. (A-D) Distribution of gene mutations differs between high- and low-risk groups in the Vizome database. (E) The different expression of HLA and immune checkpoint genes between high- and low-risk patients in training column.
TGFβ1 mRNA level in the bone marrow of patients with AML
GSEA analysis was carried out in high- and low-risk groups and the results showed that the high-risk group genetic profiles were significantly enriched in TGFβ signaling pathway (|NES|༞1, NOM p༜0.05, Fig. 5A), which was a main research topic in our laboratories currently.
To investigate the expression of TGFβ1 and TGFβR2, and examine whether the expression was correlated with disease characteristics and clinical outcomes in adult patients with de novo AML, we collected 33 AML and 15 IDA bone marrow biopsy specimens. In our study, the AML group showed elevated TGFβ1 mRNA levels compared to the IDA control group (p༜0.001), conversely, TGFβR2 mRNA levels showed lower level in AML patients (p = 0.09, Fig. 5B). Furthermore, a simple comparison of TGFβ1 and TGFβR2 mRNA levels between male and female patients showed no significant difference, neither in different age or risk stratification (Fig. 5C-E). Results of univariate and multivariate Cox proportional hazards regression analysis of TGFβ1, TGFβR2 and other clinical indicators were shown in Table 1. The high expression of TGFβR2, TGFβ1, and low levels of hemoglobin were considered as risk factors. Last, survival curves of TGFβ1 and TGFβR2 were plotted using the K-M method. Patients were divided into high- and low-level groups according to the cutoff value based on ROC curves, and we found the patients in low-level TGFβ1 and TGFβR2 group had a significantly OS than those in high-level groups (p༜0.05), respectively (Fig. 5F, G).
Exogenous TGFβ1 resists AML cells to chemotherapy
We then performed experiments to determine whether exogenous TGFβ1could enhance the malignancy of AML cells, especially drug resistance. KG-1a cells, human leukemia progenitors cell line which could response to TGFβ1, were selected as cell model for subsequent experiments. To determine the optimum concentration of exogenous TGFβ1, we set a series of concentration gradient for 48h, and determined 10ng/ml was used for subsequent experiments based on the smad2 phosphorylation levels (Fig. 6A). Simultaneously, this kind of phosphorylation could be inhibited by one of small molecular TGFβR1 inhibitors named SD208 when concentration was 1.0 um (Fig. 6B).
We next assessed how TGFβ1 affects KG-1a basic biological function, and found exogenous TGFβ1 slightly suppressed the proliferation of KG-1a and this inhibition was due to cell cycle arrest instead of apoptosis (Fig. 6C-E). Beside, the percentage of senescence cell, migration and invasion ability was increased after TGFβ1 treatment (Fig. 6F, G).
To determine whether TGFβ1 could protect AML cells from the pro-apoptotic effect of cytarabine (Ara-C) and doxorubicin (DOX), KG-1a pre-treated with or without TGFβ1 for 48h were incubated with 50um Ara-C and 3um DOX, respectively. We observed that the apoptosis of tumor cell was attenuated after additional of exogenous TGFβ1, the rate of early and late apoptosis cells were both reduced, importantly, this phenomenon could be reversed by SD208 (Fig. 6H, I).
The influence of exogenous TGFβ1 on the biological function for AML cell line. (A) Western blot analysis showing that the optimal concentration of TGFβ1 and SD-208 (B) according to the p-smad2 level. (C) Cell Counting Kit-8 (CCK8) assay (Beyotime Institute of Biotechnology, China) results indicating the impact of TGFβ1 on KG-1a proliferation at 0, 24, 48 and 72 hours. (D, E) Flow cytometry experiments showing the apoptosis and cell cycle by using the APC Annexin V Apoptosis Detection Kit (BestBio) and COULTER DNA PREP Reagents Kit (Beckman). Cell cycle phase distributions were generated by fitting DNA content histograms to applicable models using ModFit LT software for Windows (Version 5.0). (F) Senescence-associated‐β‐Galactosidase (SA‐β‐Gal) staining (Beyotime) showing the rate of senescence cell after TGFβ1 treatment for 7 days and cells that stained green‐blue were evaluated as positive senescent cells. (G) Transwell migration and matrigel invasion assays showing the increasing ability of motility and invasion. Cell suspension (1 × 105, 200 µl of serum-free medium) were seeded onto 8-mm Pore Transwell Inserts (Corning) in 24-well plates and photographed within 24h. (H,I) Flow cytometry data indicating the percentage of Annexin Ⅴ+ KG-1a cells co-cultured for 48h with Ara-C and DOX pre-treated with DMSO, TGFβ1 (10ng/ml) and TGFβ1 plus SD208 (1um).