Clinical and prognostic relevance of CXCL12 expression in acute myeloid leukemia

Patients and database

We used six publicly available datasets from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/) in this study. Of these datasets, three of them consist of expression data for bulk primary AML samples (the TCGA dataset, GSE6891 and GSE10358), two involve expression data for healthy and AML bone marrow samples (GSE30029 and GSE63270) and one includes expression data for major hematopoietic lineages during the differentiation (GSE42519), respectively. Expression data for FAB AML subtypes and patient characteristics were acquired from TCGA sequencing data. Normal mRNA expression data from Genotype-Tissue Expression (GTEx) database and pan-cancer cell line sequencing data from Broad Institute Cancer Cell Line Encyclopedia (CCLE) database were extracted through their portal websites for analysis.

Estimation of immune cell fractions

The CIBERSORT algorithm (https://cibersort.stanford.edu/index.php) was used to quantify immune cell fractions according to the expression pattern of CXCL12. More detailed information for this method is available in Newman AM et al.’s introduction (Newman et al., 2015). Immune cells calculated in this study contains B cells (including B cells memory and B cells naive), T cells (including CD8+ T cells, CD4+ naive T cells, CD4+ memory resting T cells, CD4+ memory activated T cells, follicular helper T cells, regulatory T cells and T cells gamma delta), macrophage (including M0 macrophage, M1 macrophage, and M2 macrophage), NK cells (including NK cells resting and NK cells activated), dendritic cells (including dendritic cells resting and dendritic cells activated), mast cells (including mast cells resting and mast cells activated), plasma cells, monocyte, eosinophils, and neutrophils. Multiple immune infiltration deconvolution methods including TIMER, quanTIseq, xCell, MCP-counter and EPIC were achieved from their websites.

Statistical analysis and bioinformatics

Low and high CXCL12 expressers were discriminated against based on the median expression level of the gene. Analyses of patient characteristics and survival were described previously (Xu et al., 2019). Differential expression was analyzed using ‘limma’ or ‘edgeR’ R package. The Pearson Correlation coefficient (denoted as Cor) was calculated to assess whether there was a correlation between cell infiltration and CXCL12 expression. Moreover, ‘ggplot2’ and ‘ggpubr’ R packages were applied for visualizing the results of data analysis. The single-cell correlation analysis of CXCL12 was conducted using single-cell data sequenced by Van Galen et al. (2019), and the published scRNA-seq dataset was downloaded from GEO: GSE116256. And gene ontology (GO) enrichment and Kyoto encyclopedia of genes and genomes (KEGG) analysis were conducted using the online website of STRING (http://string-db.org). One-way ANOVA was used to test for differences among at least three groups and the Wilcoxon test was used to determine differences in each two-group comparison. All statistical analyses and figure plotting were performed using the R software version 3.5.3 (https://www.r-project.org/). All differences were considered statistically significant at the level of P < 0.05.

A systematic overview of differential genes and cytokines in AML

Initially, we systematically analyzed a total of 261 cytokine and cytokine receptor genes (Supplemental Data 1) in AML patients and normal controls utilizing GSE30029, and all genes were screened in Supplemental Data 2 with the difference considered significant if adjusted p-value < 0.05 and |logFC| > 0. As obviously displayed in Fig. 1A, the CXCL12 gene was one of the most significantly different genes between AML patients and normal controls. Then, we further verified the differences in GSE63270 (Fig. 1B, Supplemental Data 3). The top 50 genes with significant differences were extracted from two analytic results respectively and then intersected. The common 11 genes included CXCL12, TNFRSF10B, IL1RAP, IL12RB1, IL3RA, CSF2RA, EPOR, TNFSF14, IL18BP, TNFRSF4 and IL2RG. Furthermore, considering the essential functions of CXCL12 in AML, we ultimately decided to study CXCL12 in detail.

CXCL12 expression is different among normal tissues and various carcinomas

To explore how CXCL12 is expressed across various tissues in healthy people, we analyzed the mRNA expression status utilizing the GTEx database (Fig. 2A, Supplemental Data 4). Apparently, CXCL12 expression was extremely low in blood and bone marrow while it was much higher among other normal tissues. It was also worth mentioning that the brain, stomach, and testis showed a relatively low level of CXCL12 expression compared with other tissues. More importantly, when analyzing different cancer cell lines (Supplemental Data 5), we found CXCL12 was lowly expressed among most hematopoietic neoplasms, especially in AML and chronic myelocytic leukemia (CML). Interestingly, cervix cancer cells showed exclusively lower CXCL12 expression levels than AML and CML, as shown in Fig. 2B.

CXCL12 expression is markedly reduced in AML and single cells

Considering our purpose, we studied the CXCL12 expression difference between healthy people and AML patients. A significant decrease of CXCL12 expression was identified in AML cells (GSE30029, n = 90) in comparison to normal cells (GSE30029, n = 31) (P = 3.2E−16, Fig. 3A). We detected a similar result in the GSE63270 cohort (n_AML = 62, n_normal = 42, P = 7.4E−06), as indicated in Fig. 3B. Additionally, we analyzed CXCL12 expression difference in eight FAB AML subtypes from M0 to M7 (Fig. S1) and CXCL12 transcription during hematopoiesis differentiation (Fig. S2), although the difference was of no statistical significance.

For further exploration, we delineated CXCL12 expression patterns in several hematological cells including hematopoietic stem cells (HSCs), multipotent progenitor cells (MPPs), granulo-monocyte progenitor cells (GMPs), common myeloid progenitor cells (CMPs), megakaryocyte-erythroid progenitor cells (MEPs), monocytes, and erythrocytes using the single-cell data sequenced by Van Galen et al. (2019), as shown in Fig. S3A. Interestingly, unlike the former result, MPPs and MEPs showed exclusively high expression of CXCL12. Also, we studied CXCL12 and CXCR4 expression status in the single-cell sequenced data and found CXCR4 expression was abundant in single cells from all patients while CXCL12 was rarely expressed (Fig. S3B).

Clinical features and molecular insights of CXCL12

The patient characteristics with respect to CXCL12 expression from TCGA were summarized (Table 1). Low CXCL12 expressed patients tended to have higher white blood cell (WBC) counts (P < 0.0001) and more percentages of blasts in BM (P < 0.001) and peripheral blood (PB) (P < 0.0001). With regard to frequently mutated genes in AML, low CXCL12 expressers were more likely to have FLT3-internal tandem duplications (FLT3-ITD) (P = 0.010) and NPM1 mutations (P = 0.015) and less likely to have mutated TP53 (P = 0.046) than those with high CXCL12 expression.

Decreased CXCL12 expression predicts adverse prognosis

Survival analyses were further performed using three independent cohorts (TCGA, GSE6891, and GSE10358), and all cohorts were divided base on optional cut-off values determined by the x-tile method. In the TCGA whole-AML cohort, CXCL12^low patients showed significantly decreased overall survival (OS) (P = 0.013, Fig. 4A) and event-free survival (EFS) (P = 0.0001, Fig. 4D). The impact of CXCL12 expression on outcome was also observed in the non-M3-AML (OS, P = 0.002, Fig. 4B; EFS, P = 0.0007, Fig. 4E) and cytogenetically normal (CN)-AML (OS, P = 0.002, Fig. 4C; EFS, P = 0.002, Fig. 4F) subgroups. This finding was further validated using GSE6891 (Fig. S4) and GSE10358 (Fig. S5) and the results were identical to the former.

The potential relationship between low CXCL12 expression and the existence of FLT3-ITD and NPM1 mutations led us to further explore the additional prognostic value of the gene in more subdivided groups. As expected, low CXCL12 expression patients showed significantly inferior survival probability in the FLT3-ITD absent (Fig. 5A, P = 0.003) and NPM1 mutated (Fig. 5B, P = 0.068, a trend for significance actually) patients, which could provide more subdivision proofs for AML based on the dichotomous stratification of CXCL12 expression.

Multivariate analysis of CXCL12 expression for OS in the TCGA cohorts

To evaluate whether there were potential interrelationships between CXCL12 expression and known clinical factors, we first conducted a univariate analysis in the TCGA dataset and then covered variables with a univariable P ≤ 0.20 in the multivariate Cox proportional hazard model. The multivariate analysis result of CXCL12 expression for OS was displayed in Table 2. Low CXCL12 expression maintained a high hazard rate for whole-AML patients (HR = 2.25, 95% CI [1.51–3.35], P < 0.0001), after adjusting age (P < 0.0001), WBC count (P < 0.0001), cytogenetic risk (P = 0.002) and mutation status of TP53 (P = 0.004), RUNX1 (P = 0.023) and DNMT3A (P = 0.012). The predictive significance also applied to CN-AML patients (HR = 2.14, 95% CI [1.15–3.97], P = 0.016), as indicated.

Profile of relationship between CXCL12 expression and immune cells infiltration

The landscape of the relationship between CXCL12 expression and immune cell infiltration was systematically evaluated by the CIBERSORT algorithm (Fig. 6). We found that B cells memory, plasma cells, T cells regulatory (Tregs), T cells gamma delta, NK cells activated, mast cells activated, and eosinophils were significantly changed between high CXCL12 and low CXCL12 groups, while B cell naive, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cell follicular helper, NK cells resting, monocyte, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, and neutrophils were not obviously altered between groups. More specifically, high CXCL12 expressed groups tended to harbor more memory B cells and plasma cells infiltration while low CXCL12 expressed groups exhibited more eosinophils infiltration. Less importantly, T cells regulatory (Tregs), T cells gamma delta, and NK cells were rarely infiltrated in both high and low CXCL12 expressed groups.

To verify the signatures of immune infiltration, we adopted multiple deconvolution methods including TIMER, quanTIseq, xCell, MCP-counter, and EPIC to cross-check the analytic results of CIBERSORT. Then, we made a Pearson correlation analysis between all infiltration results and CXCL12 expression, and totally 16 significantly correlated results were displayed in Fig. S6. Plasma cells (xCell, Cor = 0.273, P < 0.001) and B cells (quanTIseq, Cor = 0.246, P < 0.001; EPIC, Cor = 0.198, P = 0.007; MCP-counter, Cor = 0.145, P = 0.050) had a positive correlation with CXCL12 expression.

GO analysis and KEGG pathways analysis of CXCL12

GO analyses of CXCL12 for biological process (BP), cellular component (CC) and molecular function (MF) were performed on the online website of STRING (http://string-db.org) (Fig. 7), using genes differentially expressed between high and low CXCL12 expressers (|log FC| > 1.5, adjusted p-value < 0.05, Supplemental Data 6). GO term annotation showed that CXCL12 was correlated with anatomical structure morphogenesis, extracellular structure organization, animal organ morphogenesis, biological adhesion, system development, multicellular organismal process, animal organ development, cell adhesion, extracellular matrix organization, and multicellular organism development (BP); extracellular region, extracellular space, extracellular matrix, intrinsic component of plasma membrane, integral component of plasma membrane, cell periphery, plasma membrane, collagen-containing extracellular matrix, basement membrane, and collagen trimer (CC); binding, growth factor binding, signaling receptor binding, calcium ion binding, cell adhesion molecule binding, glycosaminoglycan binding, chemokine activity, heparin binding, receptor ligand activity, and receptor regulator activity (MF) (Fig. 7, Supplemental Data 7). GO analysis showed significantly high enrichment in anatomical structure morphogenesis for biological process, extracellular region for cellular component, and binding for molecular function. Moreover, CXCL12 was involved directly or indirectly in the protein digestion and absorption pathway according to the KEGG pathways analysis (Fig. 7, Supplemental Data 7).