Immune-enrichment of non-small cell lung cancer baseline biopsies for multiplex profiling define prognostic immune checkpoint combinations for patient stratification

Background Permanence of front-line management of lung cancer by immunotherapies requires predictive companion diagnostics identifying immune-checkpoints at baseline, challenged by the size and heterogeneity of biopsy specimens. Methods An innovative, tumor heterogeneity reducing, immune-enriched tissue microarray was constructed from baseline biopsies, and multiplex immunofluorescence was used to profile 25 immune-checkpoints and immune-antigens. Results Multiple immune-checkpoints were ranked, correlated with antigen presenting and cytotoxic effector lymphocyte activity, and were reduced with advancing disease. Immune-checkpoint combinations on TILs were associated with a marked survival advantage. Conserved combinations validated on more than 11,000 lung, breast, gastric and ovarian cancer patients demonstrate the feasibility of pan-cancer companion diagnostics. Conclusions In this hypothesis-generating study, deepening our understanding of immune-checkpoint biology, comprehensive protein-protein interaction and pathway mapping revealed that redundant immune-checkpoint interactors associate with positive outcomes, providing new avenues for the deciphering of molecular mechanisms behind effects of immunotherapeutic agents targeting immune-checkpoints analyzed. Electronic supplementary material The online version of this article (10.1186/s40425-019-0544-x) contains supplementary material, which is available to authorized users.


Background
Lung cancer accounts for the majority of cancer-related deaths, with almost two million diagnosed globally each year [1], and non-small cell lung carcinoma (NSCLC) representing 83% of cases [2]. Though surgical resection is the preferred treatment modality, most patients are diagnosed at advanced, unresectable stages. The TNM staging system has historically been the most widely used predictor of NSCLC survival. Adenocarcinoma (ADC) and squamous-cell carcinoma (SCC) subtypes have differing prognostic and predictive profiles [3]. As such, pathologists are mandated to distinguish subtypes, regardless of size and quality of biospecimens, ahead of targeted and personalized therapies [4]. Advances in subtyping have brought into question the requirement for TNM [5], and recent studies demonstrate that use of immunohistochemistry (IHC) cocktails and bioinformatics [6,7], provides comparable accuracy between poorly differentiated lung biopsies and large tumors [8,9].
The ability of T cells to control cancers is now widely accepted. The use of the adaptive immune system as prognostic and predictive is becoming standardized from indisputable evidence of immunosurveillance [10], and the Immunoscore (IM) outperforming TNM staging [11]. Though tumor infiltrating lymphocytes (TIL) are associated with positive outcomes, their anti-tumor activity is curbed by immune checkpoints (ICP). ICP-blockade therapies showing broad efficacy in NSCLC patients relative to standard care are now front-line treatments [12]. Differential responses to treatments has prompted rapid FDA approval of PD-L1 companion diagnostic (CDx) assays, and measures are being taken to address its heterogeneity and assay discordance [13]. From vast clinical successes from PD-1/PD-L1 targeting, numerous additional ICPs are being investigated as combinatorial targets or CDx to control cancer [14], autoimmunity [15], and numerous infectious diseases [16]. Initially categorized as exhaustion markers of functionally impaired T cells, ICPs are expressed by tumor-reactive TILs sharing tumor-antigen specificities and T cell receptor (TCR) repertoires with circulating ICP expressing T cells [17], suggesting these may identify responders to immunotherapies.
Diagnosis and staging of NSCLC is commonly established from core needle biopsy and fine-needle aspiration, however the size and heterogeneity of these specimens does not permit use of standard IM or PD-L1 assays, creating a critical need for the development of biopsy-adaptable CDx. We constructed a tissue microarray (TMA) from immune-dense regions of core needle biopsies from a baseline NSCLC cohort, and used it to profile infiltrating immune cell (IIC) subsets, ICPs, proliferation, and effector T cell markers. We find combinations that efficiently stratify patients, and validate prognostic ICP-signatures on additional cohorts. We profile ICP coexpression dynamics and ICP linkage to clinical parameters and IIC subsets, map ICP-interactors and associated pathways, and define the most prognostic combinations able to guide blockade therapies using baseline biospecimens of all sizes.

Study design
ICP were profiled using 17 lung cancer cohorts from differing origins, and using different methods: 1) at the protein expression level on a TMA created from a baseline NSCLC cohort (n = 81) (Additional file 1: Table S1; La Rabta Hospital of Tunis, Tunis, Tunisia); 2) at the whole-tumor RNA level using RNA-Seq datasets from two NSCLC cohorts from the TCGA, the LUAD (n = 504) and LUSC (n = 494) (http://www.cbioportal.org); and 3) at the whole-tumor RNA level using microarray datasets from 14 NSCLC cohorts from the GEO, the EGA and TCGA (n = 2435) Kaplan-Meier Plotter (http://kmplot.com). Additional breast (n = 5143), gastric (n = 2183), and ovarian (n = 1816) cohort datasets were from Kaplan-Meier Plotter. Written and informed consent procedures were approved by the ethics review committees and were obtained from patients prior to the collection of specimens. Clinical patient data was randomly numbered for complete anonymity. Censoring of cohort patient data was from time of diagnosis to last follow-up or death.

TMA construction
An illustration of the TMA construction is provided in Fig. 1a. Four μm cuts made using a microtome (Leica Biosystems) from all biopsies were α-CD45 stained for IHC using the Benchmark XT automated stainer with CC1 antigen retrieval buffer (Ventana Medical Systems) for 1 h. Slides were incubated with α-CD45 (1:50) at 37°C for 1 h, followed by the ultraView DAB detection kit and counterstaining with haematoxylin and bluing reagent (Ventana Medical Systems). Slides were scanned with an Olympus BX61VS microscope equipped with a VS110 slide scanner and 20x / 0.75 NA objective with a resolution of 0.3225 mm (Olympus). Images were exported and visualized using OlyVia image viewer software ver. 2.8 (Olympus) to identify CD45 + IIC-rich regions. Three to five IIC-rich regions of the biopsies were selected for 0.6 mm core transfer into the receiving TMA paraffin block using a TMArrayer (Pathology Devices). Paraffin blocks were kept at 4°C until used for TMA construction. TMA cores were press-sealed into place after incubation at 50°C for 10 min. The TMA was cooled at RT ON, and was chilled on ice ahead of being cut into 4 μm sections. Sections were floated onto 1 mm slides (Fisher Scientific), dried ON, and stored at 4°C until stained.

Multiplex-immunofluorescence
TMA sections were deparaffinized by incubation at 50°C for 1 h prior to 5 min incubations in successive baths (3x xylene, 95, 90, 70, and 50% ethanol, dH 2 O). Antigen retrieval was performed using Target Retrieval Solution, Citrate pH 6 (DAKO) as recommended by the manufacturer. Protein Block (DAKO) was applied against non-specific staining for 40 min. Slides were rinsed with PBS before incubation with primary antibody mixtures diluted in Antibody Diluent (DAKO), 0.05% Tween 20 (Fisher Scientific) ON in a humidified chamber at 4°C. Antibodies and their dilutions are in Additional file 1: Table S3. Following three 15 min PBS washes, slides were incubated with secondary antibody mixes for 1.5 h at RT (cross-adsorbed donkey α-rabbit, α-rat, or α-goat IgG (H&L) and/or goat α-mouse IgG1, IgG2a, IgG2b or IgGM specific secondary antibodies conjugated to Alexa-Fluors (405, 488, 594, 647 and 750) (Thermo-Fisher Scientific and Abcam) (1:250) Additional file 1: Table S3. Slides were washed with three 15 min incubations in PBS, and incubated in Sudan Black (1% in 70% ethanol) for 15 min. Slides were washed with dH 2 O for 5 min, and dried for 30 min before being set with ProLong gold antifade reagent (±DAPI) (ThermoFisher Scientific) under 0.17 mm coverslips (Fisher Scientific). Primary antibodies were individually detected by donkey α-host IgG (H&L) Alexa-Fluor 594 antibodies, and images were acquired using a Zeiss Axio Observer Z1 automated microscope equipped with a Plan-Apochromat 20x / 0.8 NA objective, a Zeiss HRm Axiocam and LED pulsed light illumination (Additional file 1: Figure S1d). Fluorescence minus one controls were used for potential fluorescence bleed-through between detection channels. In other control experiments, primary antibodies: 1) were not added, 2) were detected by alternative secondary antibodies, 3) were tested on a TMA containing 14 cancer cell lines (e.g., prostate, breast, ovarian, kidney, cervical cancer cells and Creation and analysis of IIC-enriched biopsy-based NSCLC TMA. a Illustration depicting TMA creation workflow. Baseline biopsies from a NSCLC patient cohort (n = 81) were paraffin embedded, and cut sections were stained using α-CD45 to demarcate IIC-dense regions then selected for TMA construction using original blocks. Cut sections from resulting TMA were then stained using MP-IF panels targeting immunerelated antigens including ICPs and IIC subsets. Slides were scanned to create super images permitting the development of algorithms computing antigens of interest and their colocalization for normalization (figure elements modified from Servier medical art). b Image representing α-CD45 IHC stained biopsies defining IIC-dense areas. c Example of MP-IF panels demonstrating α-ICP (green), α-CD3 (pink), α-CD4 (red), and α-CD8 (yellow) antibodies validated to surround DAPI-staining nuclei (blue). IIC-enriched core selection was performed by two different operators. TMA cores were randomized and TMAs were created by two operators. HRP, horseradish peroxidase; 2°ab, secondary antibody; AF, Alexa-Fluor dye; α, anti; μm, micron; mm, millimeter Jurkat), and 4) were replaced with isotype control antibodies (MOPC-31C, G155-178, MPC-11) (BD Pharminogen). MP-IF stained slides were scanned using an Olympus BX61VS microscope housing a BrightLine Sedat filter set (Semrock) optimized for DAPI, FITC, TRITC, Cy5 & Cy7, and equipped with a 20x / 0.75 NA objective with a resolution of 0.3225 mm and a VS110 slide scanner running FW-AS software (Olympus) that stitches individual images to build high resolution .vsi images.

Image analysis
High resolution images were imported into Visiomorph software (Visiopharm), where cores were identified and linked to patient numbers using an Array-Imager module. Using fluorescence intensity thresholding, algorithms were designed to define a region of interest (ROI) and calculate total core area, which was further trained to eliminate holes within tissues to correct for actual tissue-occupying areas (Additional file 1: Figure S1f ). Two independent operators used fluorescence intensity thresholding and size exclusions to create algorithms generating labels counting cells positive for biomarkers. Single marker labeling and co-labeling of dual, triple, and quadruple colocalizing markers were performed in the same way. For co-labeling, labels created for counting cells positive for multiple biomarkers were determined using the same thresholds used to identify and count single marker labeling cells. Created co-labels were also verified as accurately staining immune cells by two independent operators. Labels identifying markers were adjusted for IIC sizes, and were centered on DAPI staining when present in panels. Baseline fluorescence thresholds assigned for minimal signal to noise ratios determining positivity were used to calculate MFIs. Counts of algorithm-determined labels on cores were validated to reflect visual operator counts. Inter-rating correlations from algorithms created by independent operators was assessed to be > 75%. Each single or multi-marker label counts (e.g., totaling up to 15 marker permutations for each individual 5-color panel in case of DAPI + 4 markers) of individual cores were automated to be reassigned to patient ID numbers, and were then log-transformed and normalized to core size, prior to being merged with the clinical data for averaging of replicate core values, resulting in data from 73 patients for further analyses from .csv data file exports. High (hi) and low (lo) values were defined as being above or below mean ± SEM. Receiver operative characteristic (ROC) curves (SPSS software v.23, IBM), were used to validate that selected cutoff values corresponded to the best sensitivity and specificity any given marker. ICPs having inter-patient variability were found from a second method of analysis applied whereby values from individual cores were not averaged.

Statistical analysis
Power analysis determined that our retrospective biomarker study based on overall patient survival required a minimal sample size of n = 62 to reach a power of 0.80 at α = 0.05 (two-tailed) (G*Power ver. 3.1.9.2; Universitat Düsseldorf, Germany). Prism 6 ver. 6.01 (GraphPad) and SPSS software packages were used for statistical analysis of biomarkers with patient data. Log-rank (Mantel-Cox) tests with log-rank HR were used for K-M. A student's t test was used to compare two groups, and two-way ANOVA (with Tukey's or Bonferroni's multiple-comparisons tests) was used for multiple comparisons. Pearson correlation coefficients were calculated with two-tailed P values with 95% confidence intervals. P-values of less than 0.05 were considered to indicate a statistically significant difference. R with a collection of libraries was used for additional statistical correlation, linear regression, variance and clustering analysis, patient clinical characteristics and biomarker expression value relationships analyses. Here, expression values were log transformed towards a Gaussian distribution. Linear regression matrices were computed using the R glm function. Link functions were adapted phenotype distribution type (binomial, Gaussian, Poisson) for model compatibilities for explorations of relationships between biomarkers and clinical data. K-M calculations, cox model p-values and HR were validated using a survival model coupling survival status and months of survival post biopsy. PCA was used for coexpression analysis. Cumulative correlations for the expression of each ICP (and CD3-ICP) were calculated from their respective correlation matrices.

Prognostic signature validation and gene expression analysis
Kaplan Meier plotter was used to validate the prognostic value of the ICP signature, and to assess ICP gene expression modulation between tumors and normal tissues. Gene ID symbols were mapped to Affymetrix probes from GEO, EGA and TCGA datasets, and their mean expression was used to assess OS. For K-M, default settings were used with auto select best cutoff and best specific probes (JetSet probes). The 2017 version of Kaplan Meier plotter contains information on 54,675 genes for survival, including 2437 lung, 5143 breast, 1065 gastric, and 1816 ovarian cancer patients with mean follow-up times of 49, 69, 33, and 40 months, respectively. Multigene classifier function using default settings from KM-plotter was used to run the analysis on all ICPs simultaneously, where global ICP coexpression represents combined prognostic effects of all ICPs investigated in this study.

Protein-protein interaction network and pathway enrichment analysis
Identified biomarkers were subjected to comprehensive pathway enrichment analysis using pathDIP ver. 2.5 (http://ophid.utoronto.ca/pathDIP) (Additional files 2 and 3). Default settings were used, with extended pathway associations (combining literature curated core pathways with associations predicted using physical protein interactions with minimum confidence levels of 0.99). Lists were also used to retrieve physical protein interactions and explore biologically relevant links. IID ver. 2016-03 (http://ophid.utoronto.ca/ iid) was used to map identified biomarkers to proteins and retrieve their interacting partners. Default settings were used, and interactions among partners of query proteins, source information (detection methods, PubMed IDs, reporting databases), and tissue information (presence/absence of interactions in selected tissues) were included. Corresponding networks were visualized using NAViGaTOR ver. 3 (http://ophid.utoronto.ca/navigator) (Additional file 4). Word-cloud analysis was performed using Wordle software ver. 2014 (http://www.wordle.net).

Creation and analysis of immune cell-enriched tissue microarray
We aimed to develop a standardizable, immune-based, prognostic scoring method for biopsies. To reduce tumoral heterogeneity, a CD45-enriched TMA was constructed from baseline biopsies from a NSCLC cohort (Additional file 1: Tables S1 and S2). Figure 1a illustrates the construction of the TMA. Ahead of construction, nine random biopsy sections where stained for immunofluorescence (IF) using DAPI, α-CD45 and α-cytokeratin; verifying these for epithelial cancer and IIC-densities (Additional file 1: Figure S1a). Cut sections from all biopsies were then stained for IHC using α-CD45, defining IIC dense regions selected for TMA construction (Fig. 1b). IIC density of biopsies did not correlate with clinical parameters (P > 0.416) (Additional file 1: Figure S1b) or overall survival (OS) (P = 0.7880) (Additional file 1: Figure  S1c). All antibodies were validated independently (Additional file 1: Figure S1d and e), and TMAs were stained with five-color multiplex-IF (MP-IF) panels using a two-step, semi-automated method ( Fig. 1a and c). Algorithms calculated core areas to normalize labels identifying size-and fluorescence intensity-gated, colocalizing IICs and ICPs (Additional file 1: Figure S1f).
Cytotoxic and immune stimulation markers were investigated. Correlative studies between expression of effector markers (IFN-γ, GZMB, HLA-DR), and IIC subset infiltration of patient cores were used to demonstrated that expression of effector markers could be associated with the presence of CD8 + , CD4 + and CD20 + IICs (Fig.  2f ). IFN-γ (P = 0.0027) and HLA-DR (P = 0.0001) were positively associated with OS ( Fig. 2a and e). IFN-γ marks adaptive immune activation, and is central to anti-tumor immunity [26], and absence of HLA-DR is associated with metastasis [27]. IFN-γ localized to plasma membranes and periplasmic bursts of CD8 + TILs, and to nuclei of both TILs and epithelial cells (Additional file 1: Figure S1e), possibly explained by its rapid cellular export and nuclear localization signal [28]. GZMB and HLA-DR staining was typical, but rarely evident on TILs (Fig. 2g). HLA-DR is expressed by APCs [29], perhaps explaining it labeling cells neighboring CD8 + TILs. As prognostic factor for NSCLC, HLA-DR has been shown to identify M1 CD68 + TAMs [30]. GZMB labeled small cells, and is expressed by B cells, mast cells, keratinocytes, and basophils [31]. Altogether, these results demonstrate that proliferating Ki-67 + IIC; CD3 + , CD8 + , and CD4 + TILs; CD20 + TIL-Bs; and HLA-DR and IFN-γ are positive prognostic markers for NSCLC patients.  2 Highly proliferating, effector TIL and TIL-B densities are linked to positive prognostic of NSCLC patients. a Summarizing graph of P-values generated from K-M survival analyses of markers applied to IIC-enriched biopsy TMA, where significance indicates positive associations of IIC subsets, and proliferation and effector molecules with OS. b K-M curves (top) from Ki-67 co-labeling with CD45 + IICs or CD3 + TILs on TMA, and representative close-up IF images from cores (bottom) demonstrating co-labeling on cells. c K-M curves (top) from CD4 + and CD8 + TILs on TMA, with representative close-up IF images from cores (bottom) demonstrating their co-labeling CD3 + TILs. d K-M curves (top) from CD20 + TIL-Bs, PNAd + HEV, and CD68 + TAMs, with representative close-up IF images from cores (bottom). e Graph of average proportions of IIC subsets relative cell count (DAPI), where percentages represent IIC subset abundance relative to CD45 + IICs. Percentages are relative to CD45 content, and error bars represent mean ± s.d.. f Graph of correlations between IIC subsets and quantified effector molecules (IFN-γ, GZMB, HLA-DR). Percentages represent IIC subset attribution to effector molecule expression, as calculated from proportions of individual IIC subsets infiltrating cores expressing effector molecules. g K-M curves (top) of GZMB, IFN-γ, and HLA-DR effector markers, with representative close-up IF images from cores (bottom) of these markers and TILs. The number of patients (n) for each group is given on K-M curves, and remainder are in Additional file 1: Figure S2b. Algorithm design, normalization and analyses were performed by two independent operators. Norm., normalized; hi, high marker expression, lo, low marker expression; μm, micron; P, Log-rank test; ns, not significant; * P < 0.05; ** P < 0.01; *** P < 0.001; HR, hazard ratio (Logrank); CI, confidence interval of ratio; NA, not applicable NSCLC survival correlates with increased expression of ICP on TIL IFN-γ expression by activated TILs increases PD-L1 expression [32]. IFN-γ is also correlated with the expression of other ICPs, including BTLA [33], TIM-3 [34], LAG-3 [35], and PD-1 [36]. Since ICPs are expressed by various cell types, their usage as mono-CDx will lead to assay inconsistencies exemplified by PD-L1 [37]. Indeed, on our TMA, certain ICPs labeled numerous cells types (PD-L1, TIM-3, TIGIT, LAIR-1, CD73), whereas others almost exclusively labeled TILs (BTLA, LAG-3, PD-1, CD39, 2B4, CD57, CD26, CLTA-4) (Additional file 1: Figure S3a to e). Despite this, principal component analysis (PCA) demonstrated that relative to patients, tight clustering of ICPs and cognate CD3-ICPs indicated that they were mostly labeling TILs, and not other cells of the tumor microenvironment (Additional file 1: Figure S3f).
Global ICP expression is independent of immune density and provides a pan-cancer survival advantage In correlative analyses between global ICP or CD3-ICP expression and IIC subsets, IIC subset infiltration of patient cores were used to demonstrated that expression of ICPs and CD3-ICPs effector markers could be most associated with the presence of CD8 + , CD20 + and CD4 + IIC subsets ( Fig. 4a and b). We tested whether the IIC-density of biopsies influenced CD3 and ICP distributions. CD3 + TILs were highly correlated with CD45 + IICs (P < 0.0001, r = 0.3428), but global ICP expression was not (Fig. 4c), with the exception of CD3-PD-1, CD3-PD-L1, CD3-BTLA and CD3-LAG-3 (Additional file 1: Table S5). This also supports that ICPs are not uniquely expressed by TILs (ICP vs CD3-ICP; P < 0.001) ( Fig. 4c and Additional file 1: Figure S3a to e) [38,39]. ICPs correlating with CD3 were BTLA, LAG-3, TIM-3 and CD26, and CD73 and CD3-CD73 correlated with the ADC subtype [40] (Additional file 1: Table S5). Despite their clear effects on outcomes (Additional file 1: Figure S4), there was no correlation between treatments and ICP expression. We also observed that CD3-ICPs were inversely correlated with tumor size and extent ( Fig. 4d and Additional file 1: Table S5). K-M performed using global expression of ICP or CD3-ICP revealed that both positively correlated with OS ( Fig. 4e and f ), and global CD3-ICP expression also correlated with female gender (P = 0.0321, r = 0.0701).
To validate our findings on ICPs, we used the TCGA LUAD and LUSC RNA-Seq datasets. As observed from TMA analyses, advanced cancer patients and those deceased both had lower ICP expression ( Fig. 4g and h). Despite background noise from these whole-tumor RNA datasets, eight ADC patient ICPs were associated with positive OS (Additional file 1: Table S6). Additional cohorts from the Gene Expression Omnibus (GEO), TCGA and European Genome-phenome Archive (EGA) validated this finding for ADC patients (P = 4.4e-05) (Additional file 1: Figure S5), and grouped analyses confirmed that global ICP coexpression benefited NSCLC patients regardless of subtype (P = 1.1e-14) (Fig. 4i). Global ICP coexpression was also positively associated with OS for breast (P = 3.2e-03) and gastric (P = 1.3e-02), but not ovarian cancers (P = 1.6e-01), despite an observable trend ( Fig. 4j and l and Additional file 1: Table S7). These analyses also demonstrated a commonality of ICP expression in NSCLC and breast tumors relative to normal tissues (Additional file 1: Table S8). To validate the utility CDx profiling ICP on TILs, K-M was performed on ICP groups associated with OS or increased in expression, revealing that their prognostic value was maintained when coexpressing with CD4 or CD8 (Additional file 1: Table S9). These datasets also used to validate prognostic associations and increased expression of IIC subsets and T cell activation markers (Additional file 1: Table S10). Chromosomal locations of ICPs suggested that transcriptional regulation from common promoters is unlikely (Additional file 1: Table S11). Altogether, these results demonstrate that global ICP coexpression augments survival from different cancers, and their correlation with CD3 + TILs supports the development of multiplex CDx. Furthermore, since overall ICP expression was independent of IIC density, even patients with low infiltration may benefit from precision ICP-blockade therapies.

ICP combinations on TIL are associated with increased NSCLC survival
Using TMAs, we assessed minimal ICP combinations on TILs maximizing prognostic value (Additional file 1: Table S12). Indeed, the TIM-3/CD26/CD39 combination had a stronger association with OS than these did independently (P = 0.0139), and was superior when co-labeling with CD3 (P = 0.0051) (Fig. 5a). The positive effect on OS was maintained with ICPs and CD3-ICPs co-labeling for TIM-3/BTLA/LAG-3 combinations (P = 0.0018 to P = 0.0033), as it was for the 2B4/PD-1/CD57 combination (Fig. 5b and c). As supported by imaging (Additional file 1: Figure S6), comparisons of ICP and CD3-ICP K-M curves validated that these ICP combinations were specifically labeling TILs, and that the difference in prognostic association using duplex or triplex ICP panels was dependent on ICP combinations.
The feasibility of stratifying patients by adding individual ICP values instead of using ICP-colocalization values was also validated (e.g., TIM-3 + LAG-3, P = 0.0016; TIM-3 + BTLA, P = 0.0022; TIM-3 + BTLA+LAG-3, P = 0.0099), indicating that similar results could be attained from sequential IHC methods. However, our simplified method has less potential for antibody cross-reactions, loss of antigen and tissue integrity from harsh chemical treatments, loss of colocalization from permanent stains masking subsequent antigens, or potent spectral overlap of fluorescent signals requiring unmixing [41]. Altogether, these results demonstrate that the simultaneous detection of multiple ICPs on TILs using MP-IF panels efficiently stratifies NSCLC patients.
PCA was deployed to better define synergizing ICPs across different MP-IF panels (Fig. 6b). Proportions of variance of principle components (PC), corresponding to combined expression of each ICP group, validated that the first PC (PC1), followed by the second PC (PC2), accounted for the greatest degrees of variancerepresenting groups having differential and unrelated expression dynamics (Additional file 1: Figure S7a). K-M was computed using high vs low PC group values (Additional file 1: Figure S7b). From the TMA dataset, a group of highly expressed ICP (low PC1) was significantly associated with OS (P = 7.3 × 10 − 4 ). The relationship between PC1 and OS was increased using CD3-ICP values (P = 1.4 × 10 − 5 ). PC2 values representing the second ICP cluster did not demonstrate as clear a relationship with survival. Altogether, this analysis Fig. 6 RNA and protein conserved ICP coexpression groups ranked for NSCLC patient stratification. a-c Graphs depicting R package generated correlation studies made between all ICPs from RNA and TMA datasets to reveal ICP coexpression dynamics stratifying patients. From left to right, RNA expression of ICPs from the TCGA LUAD (n = 504) and LUSC (n = 494) patient samples (left two graph columns), were compared to that of ICP and CD3-ICP expression from all TMA dataset patient (n = 73) samples (right two graph columns). a Correlograms demonstrating ICP coexpression clustering, where black boxes demarcate most highly correlating ICP. b PCA for visualization of multi-dimensional ICP coexpression, relative to distributions patient data (blue circles), where yellow shaded PC quadrants are occupied by ICP coexpressing groups having positive associations with OS, defined by Additional file 1: Figure S7. c Mean correlations of ICP coexpression demonstrate those most abundantly expressed relative to all other ICPs in NSCLC patients. Analyses were performed using alternative software (see Online Methods) by two independent operators. PC1, principal component 1; PC2, principal component 2 revealed that the coexpressing ICP group BTLA +-LAG-3 + PD-1 + PD-L1 + most efficiently stratified patients across all datasets (Fig. 6b and Additional file 1: Table S14). The TIGIT + CTLA-4 + 2B4 + group was maintained across RNA datasets, and the TIM-3 + CD26 + CD39 + group was maintained across protein datasets.
We performed correlation analyses to determine which ICPs were most highly coexpressed. For RNA datasets, ICP ranking was TIM-3-TIGIT-CTLA-4-LAIR-1-BTLA-PD-1 (Fig. 6c). For TMA protein-derived datasets, this was BTLA-TIM-3-LAG-3-PD-1. In our comparison of four cancers, CTLA-4-TIGIT-PD-1-TIM-3-BTLA-LAG-3 were among those most increased in expression and having the greatest association with OS (Additional file 1: Tables S7 and S8). Additional file 1: Figure S8 demonstrates detection of ICPs from whole-tumor RNA to protein on TMA CD3 + TILs, where augmented ICPs may be at the forefront of the anti-cancer response, making these the best CDx and ICP-blockade targets. To determine whether coexpression dynamics could be reflected by time to effect on OS, we examined K-M curves to identify ICPs having the earliest effect on OS. For both RNA and protein datasets, ICPs with the greatest impact on OS, either alone or in combination (Figs. 3, 5 and 6), were among those having the earliest impact on OS (Additional file 1: Figure S9). Taken together, these results revealed that key ICP groups have conserved coexpression from whole-tumor RNA to protein on TILs, where discrepancies may arise from ICP expression by other cells of the tumor microenvironment also captured by whole-tumor RNA datasets. The prevailing conserved ICP subgroup (BTLA/TIM-3/LAG-3/PD-1) was most highly coexpressed and had the largest impact on OS. It is not known whether these ICPs are the first accumulating, or those persisting longest on TILs, but these are surely robust targets for combination CDx.

Redundant ICP-interacting proteins are linked with NSCLC patient survival
From the observation that ICPs positively associating with OS were increased in expression in tumor samples (Additional file 1: Table S8), we used the Integrated Interaction Database (IID) to identify 1750 key ICP-protein interactions from 40,555 possible interactions between all identified ICP-interacting proteins. Key ICP-interactors were refined for those that were 1) experimentally validated to interact with ICP, 2) redundantly interacting with more than one ICP, 3) associated with OS, and 4) had supporting evidence for their interactions in lung tissues (Additional file 1: Table S15). NAViGaTOR software was used to visualize all ICP-interactors, their characterized molecular functions, and supported interactions in lung tissues; demonstrating that 10 of the 13 signature ICP interacted with each other (Additional file 1: Figure S10, Table S16, and Additional file 4). Interaction networks were expanded to visualize defined groups from refined ICP-interactors (Fig. 7). The majority of ICP-interactors had a positive association with OS (64.6%); most of which also had increased gene expression in tumors (85.4%). The majority of ICPs in these two categories were also those ranking highest in interactions with other ICPs. Both increased in expression in tumors and associated with positive OS, BTLA and TIM-3 were observed to interact with a majority of these proteins (Fig. 7 and Additional file 1: Table S15). The pathDIP portal was used for comprehensive pathway enrichment analyses of ICP-ICP interactions and refined ICP-interactors lists ( Fig. 7 and Additional files 3 and Additional file 4), and word-cloud analysis was used to compile the most significant ICP-interactors and associated pathways (Additional file 1: Figure S11). Together, these results demonstrate that most ICP-interactors are increased in expression and are associated to positive outcome, further suggesting that ICPs are positive prognostic NSCLC biomarkers.

Discussion
ICPs were originally classified as exhaustion markers of functionally impaired T cells. Investigations of this reversible impairment have led to numerous clinical successes in cancer treatment. We were initially surprised that ICP expression on NSCLC TILs was positively associated with survival; a finding we confirmed using several additional cohorts spanning different solid cancers. When assessed in combinations, PD-1 and PD-L1 are positive prognostic markers of effector memory antigen-experienced CD8 + T cells [42]. ICP expression kinetics have been suggested to reflect CD8 + T cell differentiation kinetics rather than functional impairment [43], and as also suggested by our results, these are speculated to accumulate on TIL in an ordered fashion, led by PD-1, TIM-3, CTLA-4, LAG-3, and BTLA [44]. These represent robust CDx candidates because their prognostic/stratifying effects are also visible using whole-tumor RNA datasets. Another recent study by the Zippelius group is an additional demonstration of the rethinking of the meaning of T cell exhaustion/dysfunction in NSCLC, demonstrating that NSCLC TIL populations coexpressing several ICPs are highly clonal with a predominance of TCRs resulting from their antigen-driven expansion, that these secrete high levels of chemokines recruiting B cell and CD4+ helper cells into tumors, but most importantly, that this population is a strong predictor of robust responses to immunotherapy and overall survival [45].
We identify BTLA as the most reproducible prognostic biomarker spanning all cohorts investigated, as it: 1) predicted positive outcome from the TMA; 2) predicted positive outcome from whole-tumor RNA; 3) was most coexpressed with other ICPs across all datasets; 4) had earliest effects on OS; 5) had increased expression in tumors; 6) interacted with a majority of other ICPs and other proteins; and 7) was almost exclusively expressed by TILs. Responders to adoptive cell transfer (ACT) have increased proportions of CD8 + BTLA + TIM-3 + TIL infusion products [46], and BTLA is speculated to be the final checkpoint towards differentiation into effector T cells [47]. Accordingly, BTLA was the only ICP decreased from stimulation ahead of transfusion of autologous cultures used for successful NSCLC ACT [48,49]. BTLA may be an ideal target for ICP-blockade, because it is restricted to lymphoid tissues, and its inhibition restores TCR signaling [50]. BTLA protects TILs from apoptosis [51], and with T cell longevity estimated at over a decade [52], balanced BTLA expression may make the difference between antigen-experience and death.
Even using large biospecimens, the heterogeneity of the tumor microenvironment is the biggest challenge to finding prognostic and predictive biomarkers. We have thus developed a method for stratifying patients from limited biospecimens unsuitable for standard IM. Our restriction of analysis to immune-dense regions overcomes both size and heterogeneity of biospecimens, identifying several IIC and ICP combinations stratifying NSCLC patients. This fully-automatable combination CDx platform represents an optimal salvage method for profiling TILs from baseline biopsies ahead of personalized ICP-blockade therapies. The BTLA, TIM-3, LAG-3 and PD-1 combination on TILs was increased in expression and offered the best survival advantage. These ICPs were among those having: 1) highest correlation with any other ICP on CD3 + TILs, 2) positive association with OS at both RNA and protein levels, 3) the earliest effects on K-M curves, 4) equal impact on OS from the alternative method of analysis, and 5) decreased expression at advanced stages. These ICP may be among the first, or most persistently expressed by TILs gaining antigenic experience, as suggested by their strong correlation with TIL-Bs. This ICP subgroup represents the best CDx combination for stratifying patients using small biospecimens.
This work was in part performed to address the issues plaguing PD-L1 as CDx. Demonstrations of PD-L1 contribution to disease is challenging because it is easily inducible or constitutively expressed by many cell types. We observed that PD-L1 only stratified patients when co-labeling with CD8 or TIM-3. Likewise, despite initially described as a poor prognostic factor, PD-L1 association with TILs is linked to better outcomes in diverse cancer types [53,54], and its expression on TILs predicts response to α-PD-L1 [55,56]. Our finding that CD3-PD-L1 association with OS was affected by the alternative method of analysis confirms variability of PD-L1 expression on TIL within individual biopsies. Conversely, associations of CD8-PD-L1 and TIM-3-PD-L1 with OS was unaffected, substantiating little variability in their co-occurrences. Success of PD-L1 as CDx may thus not come down to the choice of clone, but rather from its profiling in combinations providing adequate 'immune contexture'. Like PD-L1, we find that numerous ICPs and IICs better stratify patients when profiled in combination.
Despite ICP being excellent targets for immunotherapies, they are also crucial for T cell survival. Our study does not aim to invalidate reports of ICPs as inhibitory receptors: Indeed, certain solitary ICP from whole-tumor RNA-datasets are associated with negative outcomes. Nonetheless, evidence that the majority of redundant ICP-interactors positively associate with outcomes implies ICPs have numerous important functional roles for T cells (Additional file 1: Table S17). In relation to our findings that TIL-Bs correlate with ICP coexpression and inversely correlate with metastasis, ADC clonal neoantigenenriched tumors are significantly associated to OS, have increased ICP expression, and are more sensitive to blockade therapies [57]. Specific ICP combinations may accumulate on TILs actively becoming educated against clonal neoantigens, and may protect TILs from apoptosis by slowing metabolism and differentiation kinetics. Robust MP-IF ICP CDx may identify TILs primed for tumor elimination, and the best targets for personalized immunotherapies. MP-IF ICP CDx may be also used to monitor ICP repertoires of tumor-reactive TIL expansion products for ACT. MP-IF ICP CDx created according to ICP ranking can anticipate additional ICPs arising during immunotherapies, and improve response rates to monoand combo-ICP-blockade towards their permanent adoption by mainstream oncology.

Conclusions
In this hypothesis-generating study, deepening our understanding of immune-checkpoint biology, comprehensive protein-protein interaction and pathway mapping revealed that redundant immune-checkpoint interactors associate with positive outcomes, providing new avenues for deciphering the effects of immunotherapies. We find combinations that efficiently stratify patients, and validate prognostic ICP-signatures on additional cohorts. We profile ICP coexpression dynamics and ICP linkage to clinical parameters and IIC subsets, map ICP-interactors and associated pathways, and define the most prognostic combinations that can guide blockade therapies using baseline biospecimens of all sizes.

Additional files
Additional file 1: Figure S1. Methods, antibodies and validation that IIC-density does not correlate with patient clinicopathological characteristics. Figure S2. CD45 distribution on TMA cores and K-M plots of immune cells and Ki-67. Figure S3. Five-color panels demonstrate which ICP subsets preferentially label TIL. Figure S4. Effect of applied treatments on OS of TMA patient cohort. Figure S5. Validation of effect of ICP signature on additional NSCLC cohorts and cancers. Figure S6. MP-IF ICP combination panels stratifying NSCLC patients. Figure S7. Kaplan-Meier survival analysis of principle components positively associated with OS. Figure S8. ICP coexpression ranking demonstrates ICP subset. Figure S9. Timing of effects of ICP expression on Kaplan-Meier survival curves. Figure S10. Compressed view of refined ICP-interactors presented in Fig.  7. Figure S11. Word-cloud analysis of top ICP interactors and associated pathways. Table S1. Clinicopathologic characteristics of TMA cohort. Table S2. Correlations between clinicopathological characteristics of TMA cohort. Table S3. Antibodies used in the study. Table S4. Correlations of TMA ICP MFI with immunogenicity and clinical characteristics. Table S5.
Correlations of TMA ICP counts with immunogenicity and clinical characteristics. Table S6. Primary validation of positive association of ICP expression with OS. Table S7. Secondary validation of positive association of ICP expression with OS. Table S8. ICP RNA expression in normal vs. cancer tissues. Table S9. Validation of increased positive association with OS by ICP-TIL combination. Table S10. Validation of effect of IIC expression on OS. Table S11. Chromosomal locations of profiled ICP. Table S12. Association of TMA ICP combinations from MP-IF panels with OS. Table S13. Figure 5a correlogram common ICP groupings. Table S14. Figure  S7 PC1 and PC2 groups positively associated with OS. Table S15. ICPinteractors having effects on K-M and modulated in their expression. Table S16. ICP-ICP interactors from IID. Table S17 Harvey of the CRCHUM bioinformatics department for his performing statistical validation of biomarkers, advanced bioinformatics, and assistance in interpreting TCGA and TMA datasets analyses. We thank Adolpho Faria of the CHUM Technology Transfer and Research Contracts Office, for his assembly and execution of Material transfer agreements. We recognize the generosity of the patients participating in this study, and are indebted to the Hopital La Rabta pathology department for their assistance with biopsy selection and clinical data assembly. We are grateful to The Cancer Genome Atlas (TCGA), Michael Trent Reznor, to Simon Turcotte and Jean-François Cailhier for their suggestions, and to Jacqueline Chung, Steve Trosok, and Jérôme Galon for their critical reviews. The authors declare no conflicting or competing financial interests. Availability of data and materials TCGA data and associated clinical data for all patients in this study are available at the cBioPortal for Cancer Genomics at http:// www.cbioportal.org/index.do. GEO and EGA datasets are available at http:// kmplot.com/analysis/, and database search results are available as Supplementary Data. Comprehensive pathway enrichment analysis, refined ICP-interactors from Protein-Protein Interaction analyses, and interactive NAViGaTOR networks are available as Supplementary Data. Any additional data supporting the findings of this study are available from the corresponding author upon reasonable request. The step-by-step protocols used in this study will be deposited to Protocol Exchange and be linked to the Online Methods section.
Authors' contributions AM designed and supervised the study; created the MP-IF panels and algorithms; analyzed and interpreted the data; organized the data and wrote the manuscript. ABA consented NSCLC patients for the TMA cohort, organized the clinical data, selected CD45-dense regions of biopsies, sectioned the TMA and applied MP-IF panels. DB imported all high resolution TMA images into Visiomorph, assigned cores to patient Id numbers, created algorithms, linked TMA data to clinical parameters, normalized and analyzed TMA datasets, and performed statistical analysis. LM and CC reviewed the α-CD45 IHCstained biopsies, and constructed the TMA and accompanying map. IJ performed integrative analysis, interaction network and pathway enrichment analysis and retrospective clinical validation of the effects of the ICP signature on OS of NSCLC patients. RL designed, supervised, approved and oversaw the study, and suggested analyses. IJ, A-MM-M, DB, LM, KH, and RL revised the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate Written and informed consent procedures were approved by the ethics review committees (La Rabta Hospital of Tunis, Tunis, Tunisia) and were obtained from patients prior to the collection of specimens. All experiments were performed at the Centre de recherche du Centre hospitalier de l'Université de Montréal (CRCHUM), Montréal, Canada, and the study was reviewed and approved by the local Institutional Review Board.

Consent for publication
Written informed consent was obtained from the patients for publication of their individual details and accompanying images in this manuscript. The consent form is held by the authors and is available for review by the Editorin-Chief.