IL-2 stromal signatures dissect immunotherapy response groups in non-small cell lung cancer (NSCLC)

Nanostring GeoMX Digital Spatial Profiler (DSP)

NSCLC TMA slides were stained and analysed by DSP at the Systems Biology and Data Science Group at Griffith University (Gold Coast, Australia). Morphology/visualisation markers consisted of CD3, CD68, and Pan-cytokeratin. mRNA panel consisted of 1812 curated genes from the human cancer transcriptome assay (CTA) including housekeeping and negative control probes. Slides were processed and hybridised according to the manufacturer’s instructions. ROI selection on each TMA core was performed such that 660 µm circles were segmented into cytokeratin⁺ (tumour) and cytokeratin⁻ (stroma) regions. Barcodes from these regions were collected to generate measurements per compartment. Barcodes were sequenced, mapped, and counted by NGS readout as per manufacturer’s instructions. Quality control (QC) was performed within the DSP analysis suite to remove outlying probes and collapse counts from 5 probes per gene to single measurements. This QC data was output for bioinformatic analysis.

Protein panel consisted of 68 antibodies including core panels across Human immune cell profiling, IO drug target, Immune activation, Immune cell typing, pan-tumour, cell death, PI3K/AKT panels. Slides were processed as per manufacturer’s instructions and tumour/stroma regions demarcated as above. Antibody barcodes were counted on Ncounter platform as per manufacturer’s instructions and External RNA Controls Consortium (ERCC) QC performed in DSP analysis suite prior to outputting data for bioinformatic analysis.

Multispectral Multiplex IHC

Slides were stained using validated MOTIF lung cancer PD-1/PD-L1 panel (Akoya Biosciences, US) as per manufacturer instructions on Leica Bond RX (Leica biosystems, US) at the Walter and Eliza Hall Institute (WEHI) histology core (Melbourne, Australia). Whole slide scans were performed on Vectra Polaris (Akoya Biosciences, US) by WEHI histology core, and images were spectrally unmixed in InForm (Akoya Biosciences, US) using MOTIF spectral libraries.

Image analysis

Image analysis was performed by Enable Medicine (US). Nuclear cell segmentation was performed using DeepCell, followed by segmentation dilation ^{19, 20}. Cellular protein expression levels were computed from the mean Opal fluorophore intensity for each biomarker. Cell classification was performed by gating on respective fluorescence intensity. Cells were assigned to classes: CD68⁺, (CD68⁺, PD1⁺), (CD68⁺, PD-L1⁺), CD8⁺, (CD8⁺, PD1⁺), (CD8⁺, FoxP3⁺), (CD8⁺, PD1⁺, FoxP3⁺), (FoxP3⁺, PD1⁺), PanCk⁺, (PanCk⁺, PD-L1⁺). Cell frequency was defined by counting the instances of each cell class and normalizing to the total cell count for each core. Cell interactions were defined by shared edges of the Voronoi tessellation generated from the cell centres ^{21, 22}. Cell interaction frequencies were calculated by counting the instances of an interaction and normalizing to the total number of interactions in the core. Neighbourhoods were defined using the k-nearest neighbour algorithm (KNN). For each cell, the cell types of the 10 nearest neighbours were assigned as the features of that cell. These features were then run through an unsupervised KNN algorithm and assigned to 10 clusters. The choice of 10 nearest neighbours and 10 clusters was chosen heuristically as parameters that worked well for a wide variety of datasets ^{21, 22}. For cross-sample comparisons, the frequency of each cell type/interaction/neighbourhood was normalized against the highest frequency among the cores analysed. The cross-sample comparisons were used to generate dendrograms of cell frequencies/interactions/neighbourhood. T-tests were performed between response status groups for each of these metrics.

Data analysis

Data analysis was performed by Queensland Cyber Infrastructure Foundation (QCIF, Qld, Australia). Data quality was evaluated by principal component analysis (PCA) and coefficients of variation (CV) assessed to determine suitability of the RUV-III normalisation method ^{23, 24}. Differential analysis was performed within DeSeq2 and Limma packages ^{25, 26} and sparse partial least squares-discriminant analysis (sPLS-DA) was performed within mixOmics package ²⁷. The Kaplan-Meier survival analysis (KMSA) and tree structured regression models was conducted within R studio ²⁸ using Survival, Party and Partykit package ^29–31 and plots generated by ggplot2 ³². Ingenuity Pathway Analysis (IPA®) was used to evaluate upstream regulators of differentially expressed transcripts (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis).

Patient Cohorts

The IO cohort (n=41) consisted of patients who received ICI therapy in the second line setting, of which 39% were classified as responsive (R), while 61% were non-responsive (NR) (Table 1). Anti-PD-1 therapies Nivolumab and Pembrolizumab comprised 94% of treatments, with one patient receiving anti-PD-L1 agent Durvalumab. 94% and 32% patients remained alive at follow up time, for responsive, and non-responsive groups respectively. The cohort contained squamous cell and adenocarcinoma NSCLC histology and was comprised of advanced stage III or IV disease at time of ICI treatment.

Multiplex IHC Spatial analysis

The current PD-L1 IHC companion diagnostic for PD-1/L1 axis ICI therapy remains problematic, both in its clinical assessment and biological meaning ³³. Here were utilised a robust multispectral, multiplexed IHC (mIHC) assay ³⁴ to evaluate the frequency and spatial associations of ICI targets PD-1/PD-L1, cytotoxic T-cells (CD8), T-regs (FoxP3), and macrophages (CD68) (Figure 1 A, D). Normalised cellular counts were assessed for each marker and their phenotypic subsets (e.g. CD8⁺, FoxP3⁺) (Figure 1 B, E, Supp Figure 1 A, B), their direct cell: cell interactions (Supp Figure 1 C, D) as well as nearest neighbour classifications (Figure 1 C, F, Supp Figure 1 E, F) with enrichment of each of these features in ICI response groups evaluated (Figure 1G). Neither tumour or immune cell expression of PD-L1, nor T cell expression of PD-1 was enriched in the responding cohorts, and in fact, the frequency of any phenotype alone did not associate with response in our data. Significantly however, the interactions between CD68⁺ and PD1⁺, FoxP3⁺ cells were enriched in ICI refractory tumours (p=0.012), while a trend existed for the interaction between CD8⁺, PD1⁺, FoxP3⁺ exhausted T cells and PanCk⁺, PD-L1⁺ tumour cells in ICI resistant tumours (p=0.073) (Figure 1G). K-nearest neighbour analysis (Figure 1 C, F) did not indicate enrichment of cellular neighbourhoods described by the markers evaluated in this assay (Figure 1 G, Supp Figure 1G).

• Download figure
• Open in new tab

Figure 1. Multispectral spatial analysis revealed the enrichment of cellular interactions in ICI sensitive tumours.

Representative analysis of TMA cores from responders (A-C) and non-responders (D-F). A Representative NSCLC core from responder. B Concordant cell-type Voronoi. C Concordant neighbourhood Voronoi. D Representative NSCLC core from non-responder. E Concordant cell-type Voronoi. F Concordant neighbourhood Voronoi. G Volcano plot indicating enrichment and significance of cell frequency, cell interaction and cellular neighbourhoods in responding/non-responding NSCLC cohorts.

Digital Spatial Profiling

While targeted multiplex IHC offers parallel spatial insights into tumour cellularity, it falls short of the depth required to tease apart the composition of complex tissues. Digital spatial profiling was therefore applied to better delineate the immune proteome and targeted transcriptome of the IO cohort tissues. Regions of tumour cores were segmented using cytokeratin/ non-cytokeratin immunofluorescent masks to enable transcriptome and proteomic marker measurement within tumour and stroma compartments independently (Figure 2). In this way, the measurements obtained from tumour regions comprised both tumour cells as well as infiltrating host cell populations. Similarly, stromal regions described the cellular content of tumour-associated stroma. While most cores contained both tumour and stroma regions, some punches lacked sufficient stromal content for the mRNA assay. As such, of the 41 IO cohort samples, robust transcriptome data was obtained for 33 tumour ROIs and 23 stroma ROIs, which was matched for 21 patients, while protein measurements were complete in both cohorts due to the nature of antibody-based detection compared to that of NGS. The resulting data was assessed for appropriate normalisation, where housekeeping approaches were compared to Remove Unwanted Variation (RUV-III) ²⁴ factor analysis which has been demonstrated to be more effective for variance stabilisation than traditional normalisation methods ³⁵ (Supp Fig 2). One tumour region of protein assay was excluded from analysis.

• Download figure
• Open in new tab

Figure 2. Tumour cores were segmented into PanCK⁻/⁺ regions for DSP analysis of proteomic and mRNA features.

A Representative field of non-responsive NSCLC tumour. B Representative field of responsive NSCLC tumour. C Tumour (red) and stroma (aqua) mask applied to tumour shown in A. D Tumour (red) and stroma (aqua) mask applied to tumour shown in B. Immunofluorescence markers Green = Pan-cytokeratin, Red = CD45.

Unsupervised dimension reduction was applied to explore the data’s relationship with sample clinical features, including IO response, stage at diagnosis and NSCLC histology, however variance across samples was not clearly attributable to these sample characteristics (Supp Fig 3).

Postulated markers of IO response by DSP

To further test for associations of cellular markers postulated to be involved in IO response, we evaluated therapeutic targets PD-1 and PD-L1 ^{11, 16}, tumour cell antigen presentation (HLA-DR) ^{36, 37} and cytotoxic T cells (CD8⁺) ¹² in our data. Comparative analysis by DSP alongside our mIHC data indicated that these markers were not significantly associated with ICI response in protein or mRNA DSP assays, despite a trend for higher abundance of PD-1 and CD8 mRNA in tumour regions of ICI sensitive patients (Supp Fig 4).

Differential analysis of IO response

With this in mind, we sought to identify markers in our DSP data which better discriminated ICI response. Differential analysis was performed between responders (R) and non-responders (NR) for protein markers in tumour regions (R(n)=15, NR(n)=24), stromal regions (R(n)=16, NR(n)=24), as well as mRNA in tumour (R(n)=15, NR(n)=19), and stroma (R(n)=8, NR(n)=18) regions (Figure 3, Supp Table 1).

• Download figure
• Open in new tab

Figure 3. Differential expression of protein and mRNA features by ICI response.

A Volcano plot of DE proteins in tumour. B Volcano plot of DE mRNAs in tumour. C Volcano plot of DE proteins in stroma. D Volcano plot of DE mRNAs in stroma. Tables indicating most significant DEs shown. Negative value indicates upregulation in responding patients. Differential expression was performed by IO response status within DESeq2 using RUV normalised data.

This analysis indicated the upregulation of IL-2 receptor alpha (CD25) (p=0.028) in tumour regions of responding patients (Figure 3A), which corresponded with higher expression of IL-2 mRNA (p=0.001) within their stroma (Figure 3D), suggesting key conditions for ICI efficacy. Natural killer cells (CD56) (p=0.005) and immuno-inhibitory markers CTLA4 (p=0.021), IDO1 (p=0.023) were lower in responders’ stroma (Figure 3C), with decreased levels of GITR (p=0.01), CD4 (p=0.01), CD56 (p=0.04), PD-L2 (p=0.04) as well as putative cancer stem cell marker CD44 (p=0.02) and CD24 mRNA (p=0.0008) in their tumour regions (Figure 3A, B). Interestingly, the expression of CD44 receptor in tumour regions shared an inverse relationship with the mRNA of its stromal ligand, SPP1 (osteopontin). Significantly more SPP1 mRNA (p=0.008) was observed in the stroma of responding patients (Figure 3D), while its tumour cell receptor, CD44, was relatively depleted in these patients. Decreased PTEN was associated with response in both tumour (p=0.02) and stroma (p=0.02). Surprisingly, EPCAM protein (p=0.01) was highly abundant within the stroma of responding patients, and its stromal expression appeared specific as evidenced by a lack of concordant epithelial cell marker, pan-cytokeratin.

IL2 is a key cytokine with pleotropic results on immune cell recruitment, activation, and survival. To investigate the cell phenotypes potentially being influenced by IL2 production, correlations between stromal IL2 mRNA and concordant stromal protein marker data were evaluated (Supp Fig 5). Interestingly, pro-apoptotic markers cleaved caspase 9 (R=0.73, p=2e-5) and BAD (R=0.63, p=5.5e-4) as well as the IL7 receptor (CD127) (R=0.52, p=7e-3) and IL2 receptor (CD25) (R=0.53, p=5e-3) were positively associated with levels of IL2 mRNA and thus ICI response, while levels of memory T cells (CD45RO) (R=-0.62, p=7e-4) possessed an inverse association.

ICI Survival Associations

Differentially expressed features were further examined for their prognostic roles given their associations with ICI response (Figure 4, Supp Table 2). Enrichment of tumour CD44 (HR=1.6, p=0.01) and stromal CTLA4 (HR=1.78, p=0.003) and CD56 (HR=1.58, p=0.07) markers in ICI refractory patients corresponded with poorer OS (Figure 4A). Interestingly, non-DE proteins including IL2 mRNA correlate, proapoptotic BAD (HR=0.5, p=0.01), and MDSC/M2 macrophage ARG1 (HR=2.37, p=0.01) markers were associated with improved and poorer outcome, respectively (Supp Table 2). 29 of 109 DE tumour transcripts, as well as 4 of 41 DE stromal transcripts exhibited prognostic associations (Supp Table 2). Segregation of these DE transcripts by their relative expression in ICI responsive patients indicated an inverse prognostic association (Figure 4B, C). Downregulation of several transcripts corresponded with significantly poorer outcome, including stromal E-selectin (SELE) (HR=652, p=8.8e⁻⁴) and T cell recruitment chemokine CCL17 (HR=70, p=0.006) (Figure 4B). Interestingly, while stromal MTOR expression was not associated with response, it was highly prognostic (HR=1065, p=0.008) (Supp Table 2). Conversely, enrichment of DE transcripts in ICI responsive patients corresponded with significantly enhanced OS, including stromal PLA2G2A (HR=0.03, p=0.02), and tumour NRP1 (HR=0.1, p=0.002), and NOD1 (HR=6e⁻⁴, p=0.003) (Figure 4B). Taken together, these results indicate the association of these biomarkers with both ICI response and OS outcome following treatment.

• Download figure
• Open in new tab

Figure 4. Cox proportional hazard survival analysis of differentially expressed protein and mRNA features by ICI response.

A Forest plot showing HR and 95% CI of downregulated proteins. B Forest plot showing HR and 95% CI of downregulated transcripts. C Forest plot showing HR and 95% CI of upregulated transcripts

Multivariate modelling of IO response

While differential analysis was key in resolving distinguishing features of our data relative to ICI response, their integration may yield an improved diagnostic signature. Multivariate modelling by sparse partial least squares-discriminant analysis (sPLS-DA) was thus performed to discover features that collectively discriminated patient response. sPLS-DA offers a method of feature selection to identify the most predictive or discriminative factors in a dataset that classify samples in a supervised framework, i.e. samples labelled by ICI response ²⁷. Features within protein (Figure 5) and mRNA (Figure 6) data were able to efficiently stratify ICI response within tumour and stroma compartments as demonstrated by patient sample separation on components 1 and 2 of the sPLS-DA model (Figure 5 & 6 A-B).

• Download figure
• Open in new tab

Figure 5. Multivariate sPLS-DA of protein features that distinguish NSCLC IO response.

Patient sample discrimination defined by components 1 and 2 of sPLS-DA model for tumour (A) and stroma (B). C Component 1 loading for tumour regions. D Component 2 loading for tumour regions. E Component 1 loading for stroma. F Component 2 loading for stroma. G ROC curve to evaluate component 1 tumour signature. H ROC curve to evaluate component 2 tumour signature. I ROC curve to evaluate component 1 stroma signature. J ROC curve to evaluate component 2 stroma signature. Colour of component loadings indicates patient group in which feature was maximally expressed. Positive or negative values in bar chart indicate positive or negative loading to the discriminant signature. Blue = Non-responder, Orange = Responder.

• Download figure
• Open in new tab

Figure 6. Multivariate sPLS-DA of mRNA features that distinguish NSCLC IO response.

Patient sample discrimination defined by components 1 and 2 of sPLS-DA model for tumour (A) and stroma (B). C Component 1 loading for tumour regions. D Component 2 loading for tumour regions. E Component 1 loading for stroma. F Component 2 loading for stroma. G ROC curve to evaluate component 1 tumour signature. H ROC curve to evaluate component 2 tumour signature. I ROC curve to evaluate component 1 stroma signature. J ROC curve to evaluate component 2 stroma signature. Colour of component loadings indicates patient group in which feature was maximally expressed. Positive or negative values in bar chart indicate positive or negative contribution to the discriminant signature. Blue = Non-responder, Orange = Responder.

Within tumour regions, responding patients could be characterised by increased NY-ESO-1 and lower levels of CD4, PTEN, PDL2, IDO1, CD44 and CD56 markers (AUC=0.96) (Figure 5C) with a secondary signature comprised of higher ER-alpha and lower GITR, CD25, CD14, STING and Ki-67 (AUC=0.98) (Figure 5D). Furthermore, mRNA of responding tumours were characterised by higher APLNR, PDGFC, RUNX1T1, CD33, GPX4, and lower levels of CD24, JUP, GAS1 (AUC=0.97) (Figure 6C). A secondary tumour mRNA signature comprised higher levels of KITLG, STAT3, CD47 and reduced NTRK2, PPARGC1A, CCNA1 (AUC=1) (Figure 6D).

Stromal protein markers that indicated IO response comprised higher levels of EpCAM, BAD, NY-ESO-1 and reduced CD56, CTLA4, IDO1, PTEN, ARG1, GZMA and OX40L (AUC= 0.91) (Figure 5E), with a secondary signature composed of increased ER-alpha, GITR, and decreased CD20, BIM, CD45RO (AUC= 0.981) (Figure 5F). Interestingly, IL2 RNA alone almost solely formed component 1 of the stromal signature (AUC= 0.92) (Figure 6E), while increased MMP9, PLA2G2A, IL33 and reduced COL4A6, WNT16, CRLF2, GRIA3, composed the component 2 signature (AUC=1) (Figure 6F).

Indications from the sPLS-DA together with results of the differential analysis provide evidence for a role of stromal IL2 and its cognate tumour receptor in the efficacy of PD1/L1 ICI therapy. Additionally, lower levels of natural killer cells (CD56) and IDO1 in both tumour and stroma were associated with response. Interestingly, the sPLS-DA model identified that tumour NTRK2 and stromal COL4A6 (collagen type IV) expression within the ICI resistant group contributed significantly to the discriminant signature, indicating their potential importance in refractory disease.

Network and cellular Inference

Ingenuity Pathway Analysis© (IPA®) was used to infer the molecular pathways associated with differentially expressed genes (DEGs) in responding patients. Stromal DEG’s indicated the potential suppression of IFNγ activity (activation score −2.22, p=0.0254) in responders, evidenced by the downregulation of several of its downstream transcriptional products in our data including CASP3, CCL17, DDX58, HLA-DQB1, IFITM1, PROM1, SELE, and VCAM1 (Figure 7A). Tumour DEGs also implicated the inhibition of estrogen receptor (ER) (activation score −2.891, p=0.013) and Wnt-1 (CCN5) (activation score −2.45, p=0.003) signalling pathways in responding patients (Figure 7 B-C).

• Download figure
• Open in new tab

Figure 7. IPA® assessment of differentially expressed genes to infer upstream regulatory pathways.

A IFNγ pathway activity was predicted to be inhibited within the stroma of responding patients. B Estrogen receptor signalling activity was predicted to be inhibited within tumours of responding patients. C CCN5 (Wnt-1) signalling was predicted to be inhibited within the tumours of responding patients. D Table of signalling pathways predicted to be influenced in responding patients, showing activation z-score and p-value of downstream mRNA target representation. Blue = inhibition, orange = activation, Red=mRNA upregulated, Green = mRNA downregulated

DEG enrichment within the Nanostring CTA pathway annotations showed some associations with the dysregulation of PI3K/Akt signalling, Differentiation, MAPK signalling, Cell Adhesion & Motility, Oxidative Stress, Interleukin signalling, and GPCR signalling in tumour regions (Supp Table 3). Similarly, stromal DEGs were represented in chemokine, GPCR and MAPK signalling pathways (Supp Table 3).

SPP1 mRNA (osteopontin) was the most highly enriched transcript (Log₂ FC=1.5) in ICI responsive patients. We therefore sought to assign the cellular source of this transcript in published NSCLC single cell RNAseq data that described nine distinct cell lineages in NSCLC tissue, including the immune cell subsets (T, NK, B, myeloid, MAST) in addition to epithelial cells, fibroblasts, endothelial cells (GSE131907) ³⁸. We observed overrepresentation of SPP1 within the myeloid cell lineage (cluster 1), with more specific expression within monocyte derived macrophages (cluster 5) (Supp Figure 6). The enrichment of SPP1 mRNA comes despite the absence of increased parental myeloid and macrophage cell markers in our data, including CD11c, CD163, CD68 (Supp Table 1), suggesting myeloid activation and secretion of soluble osteopontin is more directive in ICI response than cellular infiltration alone.