Expression profiling of single cells and patient cohorts identifies multiple immunosuppressive pathways and an altered NK cell phenotype in glioblastoma

Ethical approval for this study was granted by the ethics committee at the Leeds Teaching Hospitals NHS Trust, Leeds UK (REC number 10‐H1306‐7).

CIC classificationwas applied to GBM tumour transcriptome data from The Cancer Genome Atlas (TCGA). Briefly, consensus cluster analysis of melanomas used the expression of 380 genes specific to 24 immune cell types (; this produced six subtypes which we termed CICs. The average expression of each gene within each CIC is the cluster centroid. Using TCGA data, we used the nearest centroid methodto classify each GBM tumour into one of the CICs according to the highest Spearman’s correlation coefficient with the centroids. For each of the 24 immune cell types of the immunome compendium, we calculated a score per GBM tumour, graphically represented using a heatmap. 14 15 16 15

To compare the expression of selected genes in different patient groups we used RNAseq data from TCGA [obtained via The Cancer Imaging Archive (TCIA) 7], assigning patients to either CIC2 or CIC4 (as above). In addition, we used microarray data from the REMBRANDT study 17 downloaded from Betastasis.com. For REMBRANDT, patient samples were classified as granzyme A (GZMA)^high expressors or GZMA^low expressors based on the median expression value (n = 214). Expression of selected genes was compared between CIC2 and CIC4 (for TCGA data) or the GZMA^high and GZMA^low patients (for REMBRANDT) and analysed using non‐parametric, unpaired statistical testing (using GraphPad Prism).

Single‐cell (sc)RNAseq data were downloaded from 18 and expression of candidate genes analysed. Data were visualized using rStudio version 1.0.143 (package: gplots 3.0.1) using heatmap.2. Untransformed data clustering (unsupervised) was performed (Euclidean distance). Individual cells were classified as ‘tumour’ or ‘immune’ according to co‐expression of SRY‐box transcription factor 9 (SOX9) and epidermal growth factor receptor (EGFR) 18 and protein tyrosine phosphatase receptor type C (PTPRC), respectively. For all genes, expression of > 0 was scored positive. For immune (PTPRC⁺) and tumour (SOX9⁺EGFR⁺) cells, the number of different immunomodulatory molecules expressed was counted for each cell and the percentage of immune and non‐immune cells expressing immunomodulatory genes plotted. Expression of individual genes in non‐tumour‐bearing brain tissue was downloaded from 19. These data are also available (with graphical output) at BrainRNAseq.org.

After ethical approval and informed consent, tumours were resected and stored in phosphate‐buffered saline (PBS) or within the cavitron ultrasonic surgical aspirator (CUSA) 20. Samples were washed in PBS, CUSA samples were prepared as shown by Schroeteler et al. 20 and all samples were filtered through a 40‐µm cell strainer, washed twice in PBS, centrifuged at 400 g for 5 min and resuspended in PBS, 0·5% bovine serum albumin (BSA) and 0·05% sodium azide. Matched patient blood was diluted with PBS, layered over Ficoll (Axis‐Shield PoC, Oslo, Norway) and centrifuged at 800 g for 20 min. Tumour and blood‐derived cells were stained with appropriate antibodies and isotype controls (see Supporting information, Table S1), with single stain controls on tumour samples used for compensation during analysis using the cytexpert compensation matrix. All samples were run on a CytoFlex S (Beckman Coulter Life Sciences, Indianapolis, IN, USA) (see Supporting information, Table S1). Gated, isotype control stained, intratumoral or peripheral blood NK cells from each patient (Supporting information, Fig. S1) were assigned a gate of 2% positive, and specific antibody staining is reported within this gate.

Neural progenitor cells (NP1) were isolated from a patient undergoing surgery to treat epilepsy. The primary lines, GBM1 and NP1, were generated at the Scripps Institute. GBM11, GBM13 and GBM20 were derived at the University of Leeds using the same method and culture conditions. Peripheral blood mononuclear cells (PBMC) were isolated from whole blood of healthy donors as above. NK cells were further separated using an NK cell isolation kit (Miltenyi Biotec, Bergisch Gladbach, Germany), and cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 10% human AB serum (Sigma‐Aldrich, Gillingham, UK). 21 22

GBM stem‐like cell (GSC) lines were harvested using 0·25% trypsin/ethylenediamine tetraacetic acid (EDTA) and fluorescently labelled for 60 min at 37°C and 5% CO₂ in serum‐free media with one of the the following cell dyes: 0·4 μM cell tracker^TM (CT)‐green CMFDA (488 nm excitation), 2 μM CTorange‐CMRA (488 nm excitation), 2 μM CTviolet‐BMQC (407 nm excitation) or 5 μM calcein blue‐AM (407 nm excitation) (all from Invitrogen, Carlsbad, CA, USA) All populations were washed three times, mixed together and plated at a density of 1 × 10⁶ total cells/well in 96‐well round‐bottomed plates (Nunc, Roskilde, Denmark). Cells were stained as per the manufacturer’s instructions with 242 antibodies from the BD Bioscience Lyoplate screening panel, followed by Zombie NIR (Biolegend, San Diego, CA, USA) for 30 min before resuspension and analysis by flow cytometry. Cells were gated based on their emitting fluorescence at 520 nm (CTgreen loaded), 580 nm (CTorange loaded), 540 nm (CTviolet loaded) or 449 nm (calcien blue loaded). The median fluorescence intensity (MFI) for each gated population, for each antigen and isotype control emission at 668 nm (Alexa647 emission) was generated and GSC lines scored as positive if more than 20% of the population expressed the antigen. Flow cytometer and settings are as described earlier; analysis was performed using FacsDiva (BD Biosciences, San Jose, CA, USA), FlowJo (Treestar, Inc., Ashland, OR, USA) and Kaluza (Beckman Coulter) software.

Target tumour cell lines were labelled with the relevant cell dye (see surface screen) for 1 h at 37°C, washed twice and plated at 2 × 10⁵/well. NK cells were pre‐activated with 20 ng/ml interleukin (IL)‐15 for 48 h and mixed with targets at the E : T ratios indicated. After 5 h, cells were pelleted (300 g for 5 min), washed with PBS and stained with Zombie NIR (Biolegend) for 15 min at room temperature. Competitive cytotoxicity assays were set up as above; the two target cell types under test (GBM and neural progenitors) were labelled with either CTgreen or CTviolet, mixed 1 : 1 and used as a target population at an E : T of 5 : 1.

Effective therapy for GBM will require the elimination of the radioresistant GSCs that are largely responsible for recurrence. While tumour‐associated antigen‐specific T cells offer a highly selective therapeutic approach, antigen‐independent effector cells, such as NK cells, have the potential to target and destroy GBM tumour cells that have a low neoantigen load. 23

We used three patient‐derived GSC lines 24 shown to exhibit a stem cell‐like expression profile and recapitulate high‐grade gliomas in orthotopic xenograft mouse models 22, 24 and performed cytotoxicity assays using peripheral blood‐derived, IL‐15‐activated NK cells to confirm NK cell‐mediated killing. Tumour cells differentiated from GSCs are more sensitive to NK cell killing than the GSC themselves 26, but GSCs are killed by NK cells in the presence of activating cytokines (Fig. 1a) 25. We further tested whether NK cells activated with IL‐15 would be efficient killers of GSCs, but retain specificity for GSCs over normal neural progenitor cells. We performed an NK cytotoxicity assay using a mixed target cell population comprised of tumour GSC cells and normal neural progenitor (NP) cells 22 at a ratio of 1 : 1. For all donors, IL‐15‐activated NK cells killed tumour cells in preference to the NP cells (Fig. 1b). These results suggest that, in short‐term in‐vitro cultures when sufficient immune cells are present and activated, GSCs are an effective and preferential target for NK cells.

The presence of infiltrating NK cells in GBM 11, coupled with their ability to recognize and kill GSCs (Fig. 1), suggests that they are rendered non‐functional in the GBM tumour microenvironment. We performed flow cytometry‐based analysis of intratumoural NK from GBM tissue and compared their surface phenotype to NK cells derived from autologous peripheral blood as well as blood from healthy donors. NK cell populations were defined as NKp46⁺ and CD3^– due to high expression of CD56 (NCAM‐1) on GBM tumour cells within the sample (Supporting information, Fig. S1). To confirm sampling of immune cells from within the GBM tumour tissue (and not from blood contamination of the tumour sample) we assayed the expression of PD‐1 on T cells, and showed significantly enhanced expression of PD1 on tumour‐derived T cells compared to their blood counterparts (Fig. 2a). Expression levels of NK cell surface molecules were similar on the blood‐derived NK cells from both healthy donors and GBM patients. However, the expression of the tumour‐sensing NK cell‐activating receptors NKp30, NKG2D and DNAX accessory molecule‐1 (DNAM‐1) and the surface molecules tetherin/CD317 and CD2 were all significantly reduced on the GBM tumour‐derived NK cells compared to those from matched peripheral blood (Fig. 2b). Together with higher expression of PD‐1 on GBM‐derived T cells compared to matched peripheral blood (Fig. 2a), we also found higher expression of lymphocyte‐activation gene (LAG‐3) and CTLA‐4 (although differences in CTLA‐4 expression did not reach statistical significance) (Supporting information, Fig. S2).

Members of the TGF‐β family are highly expressed in GBM, and are important in maintaining the GSC pool 25. Furthermore, we and others have previously shown that TGF‐β reduces the expression of NKp30, NKG2D and DNAM‐1 on NK cells and is associated with their functional inactivation 26, 27. Importantly, TGF‐β induces the expression of the tetraspanin CD9 on NK cells 28, and we detected significantly increased expression of CD9 on the surface of the GBM‐resident NK cells compared to NK cells from matched peripheral blood (Fig. 2c). The reduced expression of NK cell activating receptors coupled with the increased expression of CD9 is suggestive of TGF‐β‐mediated evasion of NK cell cytotoxicity in the GBM microenvironment.

Collectively, GBM resident immune effector cells clearly demonstrate two separate phenotypes: the reduced expression of NK cell activation receptors and the increased expression of immune checkpoint molecules on T cells.

The GSC lines are selectively targeted by NK cells in vitro, but evade NK cells and other immune effector cells in vivo. To understand which immunomodulatory molecules expressed by GSC might be responsible for immune activation and inhibition, we analysed GSCs for the expression of cell surface immunomodulatory molecules. Using a flow cytometry‐based screen, we identified 116 cell surface antigens expressed on four GSC lines lines (Supporting information, Table S1). Molecules detected on the GSCs included those associated with the cancer stem cell phenotype (CD24, CD44 and CD90) (Fig. 3a), as well as widely expressed cell surface molecules, such as major histocompatibility complex (MHC) class I (and β2‐microglobulin), CD71 and CD98, as expected. Several immune inhibitory molecules were highly expressed, such as the immune checkpoint ligands programmed cell death ligand 1 (PD‐L1) (CD274) and PD‐L2 (CD273), providing a source for inhibition of PD‐1 expressing T cells (Fig. 3a). In addition, we found expression of the ectonucleotidase CD73 that, together with CD39, generates extracellular adenosine to inhibit both NK cells and T cells via purinergic receptors 29, as well as expression of CD200 and CD47, modulators of myeloid cell activity (Fig. 3a). Ligands of NK cell activation receptors, such as MICA/B (NKG2D ligand) and CD112 (a DNAM‐1 ligand), as well as CD80 (a T cell co‐stimulator), were detected, together with CD54 [intercellular adhesion molecule (ICAM)‐1] and CD50 (ICAM‐3); ligands of lymphocyte function‐associated antigen 1 (LFA‐1) required for NK cell and T cell‐mediated cytotoxicity 30. The GSC cell surface screen therefore revealed expression of a repertoire of targetable cell surface molecules with the potential to activate and inhibit NK cells, T cells and myeloid cells. This prompted us to explore the expression of immunosuppressive pathways in more detail, using a publicly available GBM single‐cell gene expression data set 18. Among 3589 single cells, we identified 757 co‐expressing SOX2 and EGFR (defined by Darmanis et al.18 as tumour cells) and 1527 cells expressing PTPRC (encoding CD45, a marker of cells of haematopoietic origin.) We next performed unsupervised hierarchial clustering using expression of lineage marker genes and genes encoding candidate immunosuppressive functions, which identified two main groups: non‐immune (comprising tumour and neuronal cells) and immune cells (PTPRC⁺) (Fig. 3b). The immune cell group was dominated by expression of numerous myeloid cell markers (Fig. 3b). Genes encoding immunosuppressive functions were expressed within both the immune and non‐immune clusters (Fig. 3b) and, overall, individual immune cells expressed a greater number of immunosuppressive genes than tumour cells (Supporting information, Fig. S3A). Consistent with the altered cell surface phenotype of GBM‐resident NK cells (Fig. 2a), we found widespread expression of TGFB family transcripts accounted for by TGFB1 expression in the myeloid cells and TGFB2 and TGFB3 expression in non‐immune cells. Furthermore, human leucocyte antigen G (HLA‐G) (which plays a key role in regulating NK cell activity in pregnancy and cancer 31, 32) was also widely expressed. The HLA‐G protein inhibits myeloid cells via receptors leucocyte immunoglobulin‐like receptor subfamily B member 1 (LILRB)1 and LILRB2 31, both of which were expressed in the immune compartment at the mRNA level. We identified strong expression of the receptor‐ligand pair hepatitis A virus cellular receptor 1 (HAVCR1) (TIM3) and lectin, galactose binding, soluble 9 (LGALS9) in the myeloid cluster (92% of immune cells expressed HAVCR2 or LGALS9 and 60% expressed both genes; Supporting information, Fig. S3B). This scRNA‐seq data along with the GSC surface antigen screen shows that both tumour and immune infiltrating cells express receptors and ligands that together constitute a complex network of immunosuppression. For example, the combined action of CD73 and CD39 generate immunosuppressive adenosine 29; our data show expression of CD73 by the tumour cells (Fig. 3a) and 5'‐nucleotidase ecto (NT5E) (encoding CD39) by the immune fraction (Fig. 3b), with Mohme et al. demonstrating CD39 expression by GBM‐infiltrating T cells 13.

Furthermore, to assess whether this immunosuppressive network was induced in response to tumour, we analysed gene expression data derived from normal brain tissue(Supporting information, Fig.). Microglia/macrophages (the only cell population in the normal brain expressing PTPRC/CD45) constitutively express several immunosuppressive genes, including immune checkpoints V‐set immunoregulatory receptor (VSIR) and HAVCR2 and checkpoint ligands LGALS9, CD274 and PDCD1LG2. Some components of the immunosuppressive network found in brain tumours are therefore present in the healthy brain. 33 S4

We have previously used tumour transcriptome data to cluster melanoma patients according to their immune cell infiltrate 14. This approach identified six CICs, with one cluster enriched in cytotoxic cells (CIC2) and another (CIC4) having low immune infiltrates and significantly worse survival 14. We used this approach to classify GBM transcriptome data (from TCGA) and, like the situation in melanoma, the cohort of 154 patients clustered into the six CICs (Fig. 4a), with two main clusters CIC2 (high immune infiltrate) and CIC4 (low immune infiltrate). We found that CIC2 was significantly enriched for tumours of the mesenchymal subtype 34 (Supporting information, Table S2) that has been previously shown to have prolonged survival 21. However, unlike melanoma 14, immune infiltration (reflected in the CIC clusters) was not associated with significant differences in survival in GBM (Supporting information, Fig. S5). There was also no significant difference in mutation burden, a surrogate of neoantigen load reflecting immunogenicity 8, 9, between CIC2 and CIC4 in GBM (Fig. 4b). These data demonstrate that patients can be stratified based on the immune infiltrate but that, unlike melanoma, this stratification has no effect on patient outcomes under the conditions of treatment currently employed.

Immune activation induces feedback inhibitory pathways, including the expression of immune checkpoint molecules, and we therefore attempted to use the expression of genes in these pathways to understand the immune environment within the GBM CIC clusters. To do this we compared expression of anti‐tumour effector functions and immunomodulatory genes in CIC2 and CIC4. Expression of the granzyme B (GZMB) and interferon (IFN)‐γ genes were significantly higher in CIC2 than CIC4, consistent with the increased infiltration of cytotoxic T cells and NK cells (Fig. 4c). Furthermore, genes encoding immune checkpoint molecules (CTLA‐4, PDCD1, HAVCR2, BTLA and VSIR), their ligands (PDCD1LG2, LGALS9) and forkhead box protein 3 (FoxP3) were also expressed at significantly higher levels in CIC2 than CIC4 (Fig. 4d), as were genes encoding soluble mediators of immunosuppression such IL‐10, TGF‐β1 and IDO1 (Fig. 4e). To confirm this we used microarray data from the REMBRANDT study 17 and GZMA gene expression as a simple surrogate for immune infiltration 35. This analysis confirmed the significantly higher expression of multiple immunosuppressive functions in patients with increased expression of anti‐tumour effector functions (Supporting information, Fig. S6). Collectively, these data drive our understanding of the GBM immune microenvironment, demonstrating a spectrum of immune infiltration, functionally compromised by an active immune‐inhibitory network.