Surface CD52, CD84, and PTGER2 mark mature PMN-MDSCs from cancer patients and G-CSF-treated donors

Abstract

Precise molecular characterization of circulating polymorphonuclear myeloid-derived suppressor cells (PMN-MDSCs) is hampered by their mixed composition of mature and immature cells and lack of specific markers. Here, we focus on mature CD66b⁺CD10⁺CD16⁺CD11b⁺ PMN-MDSCs (mPMN-MDSCs) from either cancer patients or healthy donors receiving G-CSF for stem cell mobilization (GDs). By RNA sequencing (RNA-seq) experiments, we report the identification of a distinct gene signature shared by the different mPMN-MDSC populations under investigation, also validated in mPMN-MDSCs from GDs and tumor-associated neutrophils (TANs) by single-cell RNA-seq (scRNA-seq) experiments. Analysis of such a gene signature uncovers a specific transcriptional program associated with mPMN-MDSC differentiation and allows us to identify that, in patients with either solid or hematologic tumors and in GDs, CD52, CD84, and prostaglandin E receptor 2 (PTGER2) represent potential mPMN-MDSC-associated markers. Altogether, our findings indicate that mature PMN-MDSCs distinctively undergo specific reprogramming during differentiation and lay the groundwork for selective immunomonitoring, and eventually targeting, of mature PMN-MDSCs.

Keywords: G-CSF; PMN-MDSCs; biomarkers; cancer patients; neutrophils.

Graphical abstract

Figure 1

mNDNs and mLDNs from GDs display transcriptomic similarities to mPMN-MDSCs from NSCLC and HNC patients Circulating mLDNs and mNDNs from GDs (n = 8–11), mPMN-MDSCs and autologous mNDNs from NSCLC (n = 3) and HNC patients (n = 5), as well as group-related control mNDNs from HDs (n = 10 for GDs; n = 3 and n = 4 for NSCLC and HNC patients, respectively) were isolated and subjected to RNA-seq experiments. (A) Dendrogram of the unsupervised hierarchical clustering analysis based on the top 500 highly variable genes among the GD, NSCLC, and HNC cellular sample cohorts. (B) PCA scatterplot of the top 500 highly variable genes. (C–E) Heatmaps displaying the expression patterns of DEGs identified by comparing GD (C), NSCLC (D), and HNC (E) mPMN-MDSC sample cohorts and group-related mNDNs from HDs. Data are shown as Z-scaled fragment per kilobase of transcript per million mapped reads (FPKM) values (Z score), while the clusters obtained by K-means clustering are reported on the left of each heatmap. (F and G) Heatmaps displaying the Jaccard similarity index, illustrating the extent of the overlap among the DEG clusters identified in mNDNs/mLDNs from GDs (GD cohort) and those of mPMN-MDSCs from NSCLC (F) or HNC (G) patients (NSCLC and HNC cohort, respectively). Box colors correlate with the Jaccard score value according to the Jaccard index scale, while box numbers indicate the p value for each comparison according to Fisher’s exact test. See also Figure S2 and Data S1.

Figure 2

Identification and characterization of a gene signature and associated TF regulons and upstream regulators shared by mPMN-MDSCs from NSCLC and HNC patients and mNDNs/mLDNs from GDs (A) Venn diagram showing the genes shared between cluster 1 (c1) and c2 of the GD cellular sample cohort and c1 of NSCLC and/or HNC cellular sample cohorts. Intersections highlighted in pink were selected for the generation of the mPMN-MDSC gene signature. (B) Heatmaps displaying the expression of the genes included in the mPMN-MDSC gene signature in either mNDNs/mLDNs from GDs (left), or mPMN-MDSCs from NSCLC (center) or HNC (right) patients. Data are represented as log₂ (fold change) relative to the mNDNs from HDs. Selected genes are listed on the right of the heatmap; font colors mean that they are included in the following biological process: black, metabolism; green, arginine and proline metabolism; blue, ATP metabolism; gray, hypoxia; purple, glycolysis; orange, pro-angiogenic; light blue, EMT; red, PMN-MDSC-associated genes from the literature. (C) Graph depicting the GO terms significantly (false discovery rate [FDR] < 0.05) over-represented in the mPMN-MDSC gene signature. Bars indicate the number of genes composing the top 15 enriched GO terms. (D) Heatmap displaying the activity score of the most consistently active TF regulons across all replicates of our GD and cancer cellular samples, as inferred using DoRothEA by applying the statistical framework VIPER to the mPMN-MDSC gene signature. (E) Graph showing the most significant upstream regulators of mPMN-MDSC gene signature expression in mPMN-MDSCs from HNC and NSCLC patients and GDs, as revealed by IPA. The bars indicate the p value on a logarithmic scale according to Fisher’s exact test. See also Figure S3 and Data S2.

Figure 3

The mPMN-MDSC gene signature is enriched in published signatures obtained by bulk RNA-seq of various PMN-MDSCs and by scRNA-seq of TANs from NSCLC patients and mNDNs/mLDNs from GDs (A) Enrichment of the mPMN-MDSC gene signature in published gene signatures of PMN-MDSCs^, and mLDGs from SLE (SLE mLDGs⁹) and suppressive CD16^highCD62L^low neutrophils (CD62L^low PMNs) from LPS-treated donors by bulk RNA-seq. Box colors indicate the p value calculated by Fisher’s exact test according to the reported scale: white, non-significant enrichment (p > 0.05); shades of blue, significant enrichment (p < 0.05). (B and C) UMAP representations of the integration of the TAN and NAN datasets from Salcher et al. or the TAN dataset from Zilionis et al. For the Zilionis et al. dataset, the TAN clusters defined in the original study (hN1-5) are reported. For every cell, the NAN signature score (B, left) and the TAN signature score (B, right) were calculated by evaluating the expression levels of genes previously found to be highly specific to a NAN or TAN phenotype (i.e., Figure 6H in Salcher et al.²⁹). By performing the analysis described in (B), we were able to define two main neutrophil populations defined in this study as NAN-like cells (green) and TAN-like cells (violet) (C) (see also Figures S3B and S3C). (D) UMAP plot of the mPMN-MDSC gene signature score, calculated as normalized average expression of the genes composing the signature in each cell. (E) Violin plot displaying the mPMN-MDSC gene signature score across the two subsets: NAN-like and TAN-like cells. (F) UMAP plot of the neutrophils with a high maturity score (defined as described in Figures S3E and S3F), either merged (top left) or divided by cellular sample cohort: mNDNs from HDs (top right), mNDNs from GDs (bottom left), and mLDNs from GDs (bottom right panel). (G) UMAP plot of the mPMN-MDSC gene signature score, calculated as described in (D). (H) Violin plot displaying the mPMN-MDSC gene signature score across the different cellular sample cohorts. ∗∗∗∗p < 0.0001 by Kruskal-Wallis test followed by pairwise Wilcoxon test. See also Figure S3 and Data S3.

Figure 4

Comparison of mPMN-MDSC gene signature expression between neutrophil precursors from the peripheral blood of GDs and those from BM aspirates of HDs RNA-seq experiments with neutrophils at different stages of maturation isolated from BM aspirates from HDs (PM, promyelocyte; MY, myelocyte; MM, metamyelocyte; BC, band cell; SC, segmented cell; n = 3), HD blood (mNDNs, n = 10), or GD blood (PMs/MYs, MMs, BCs, n = 4; mLDNs, n = 8; mNDNs, n = 11) were performed. (A) PCA plot based on the top 500 variable genes in mature neutrophils and their precursors purified from the peripheral blood of GDs (blue gradations) or the BM (pink gradations) or the peripheral blood (gray) of HDs. Arrows represent the putative maturation trajectories of HD or GD neutrophils. (B) Heatmap displaying the mRNA expression of the mPMN-MDSC gene signature in neutrophil precursors and mature neutrophils purified from the BM and the peripheral blood of HDs or GDs. Data are shown as Z score (FPKM). The major gene groups (g1–g3) identified by explorative K-means clustering analysis are shown. (C) Violin plots displaying the score of g1 genes by scRNA-seq data across neutrophils at different maturation stages of the GD (iLDNs, mLDNs, and mNDNs, left) and NSCLC patient (iPMN-MDSCs, mPMN-MDSCs, and mNDNs, right) cohorts. scRNA-seq data of control HD neutrophils at different maturation stages (BM-iLDNs and HD-mNDNs) are also included. The score of g1 genes was determined, in each cell, by calculating the normalized average expression of the genes that compose the g1 group. ∗∗∗∗p < 0.0001, by Kruskal-Wallis test followed by pairwise Wilcoxon test. (D–F) Graphs depicting the top 10 GO terms (FDR < 0.05) significantly enriched in g1 (D), g2 (E), and g3 (F). Bars indicate the number of genes composing the enriched GO terms. See also Figures S2 and S4.

Figure 5

Expression levels of CD52, CD84, PTGER2, ORL1, and LAIR1 mRNAs in GD, NSCLC, and HNC mPMN-MDSC samples Expression levels of CD52 (A), CD84 (B), PTGER2 (C), OLR1 (D), and LAIR1 (E) mRNAs in GD (left column of graphs), NSCLC (center column of graphs), and HNC (right column of graphs) cellular samples (n = 3–6). Gene expression is depicted as mean normalized expression (MNE) units after normalization to PPIB mRNA (mean ± SEM). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗∗p < 0.0001 by ordinary one-way ANOVA test or Mann-Whitney test.

Figure 6

Expression of surface CD52, CD84, and PTGER2 on mPMN-MDSCs and iPMN-MDSC populations from GDs and cancer patients and on mature and immature LDNs from the BM of HDs Surface CD52 (A, D, G, and J), CD84 (B, E, H, and K), and PTGER2 (C, F, I, and L) expression levels were evaluated by flow cytometry on mNDNs or mLDNs/iLDNs from GDs (n = 20) and group-related mNDNs from HDs (n = 17–18) (A–C), mPMN-MDSCs/iPMN-MDSCs and autologous mNDNs from cohorts of patients with NSCLC (D–F, n = 13–19) or HNC (G–I, n = 19–25), and mLDNs and iLDNs from BM aspirates from HDs (n = 6–8) or control peripheral blood mNDNs from HDs (HD-mNDNs, n = 5–6) (J–L). Graphs show CD52 (A, D, G, and J), CD84 (B, E, H, and K), and PTGER2 (C, F, I, and L) Δ median fluorescence intensity (MFI) values calculated as described in STAR Methods. CD52 (A), CD84 (B), and PTGER2 (C) expression on reference mNDNs from HDs is reported only for the GD cohort because representative of reference mNDNs from HDs from all other cohorts. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001 by Kruskal-Wallis or Friedman test followed by Dunn’s multiple comparison test. See also Figures S5 and S6.

Figure 7

Expression of surface LOX-1 and LAIR-1 on mPMN-MDSCs and iPMN-MDSC populations from GDs and cancer patients and on mature and immature LDNs from the BM of HDs Surface LOX-1 (A, C, E, and G) and LAIR-1 (B, D, F, and H) expression levels were evaluated by flow cytometry on mNDNs or mLDNs/iLDNs from GDs (n = 16–17) and group-related mNDNs from HDs (n = 17) (A and B), mPMN-MDSCs/iPMN-MDSCs and autologous mNDNs from cohorts of patients with NSCLC (C and D, n = 24–27) and HNC (E and F, n = 13–35), and mLDNs and iLDNs from BM aspirates from HDs (n = 6–8) or control peripheral blood mNDNs from HDs (HD-mNDNs, n = 5–7) (G and H). Graphs show LOX-1 or LAIR-1 Δ MFI values. Each symbol stands for a single donor sample. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001 by Kruskal-Wallis or Friedman test followed by Dunn’s multiple-comparisons test.