Multiplex Immunohistochemistry and High-Throughput Image Analysis for Evaluation of Spatial Tumor Immune Cell Markers in Human Breast Cancer: Preliminary Results from the Nashville Breast Health Study

Purpose Tumor-inltrating lymphocytes (TILs) have emerged as predictive biomarkers for cancer prognosis. The clinicopathological signicance of spatial TIL subpopulations is not well studied due to lack of high-throughput scalable methodology for large population studies. Methods We established a uorescent multiplex immunohistochemistry (mIHC) method coupled with computer-assisted high-throughput quantitative analysis. We then evaluated associations of six TILs markers (CD3, CD8, CD20, CD56, FOXP3, and PD-L1) with breast cancer prognosis among 188 tumor samples from the Nashville Breast Health Study. Results Our 5-plex uorescent mIHC workow was reliable, highly sensitive, non-interfering, and balanced labeling for three biomarkers per tissue section, which is applicable for high-throughput imaging quantication of spatial TILs in regular laboratory settings. In this study we found that the increased intratumoral CD56+ cells were associated with worse overall survival (OS) (HR, 4.89; 95% CI: 1.26-18.95, highest vs lowest tertiles; P trend = 0.022), suggesting the subset of immunosuppressive NK cell phenotypes within tumor bed. Increased stromal PD-L1+ cells (HR, 0.01; 95% CI, 0.00-0.89; P trend = 0.042) and a high stromal CD8+/FOXP3+ ratio (HR, 0.00; 95% CI, 0.00-0.12; P trend = 0.036) were associated with more favorable OS in stage III-IV breast cancer patients. Conclusion We established a reliable 5-plex uorescent mIHC and showed that CD56+, PD-L1, and CD8+/FOXP3+ ratio may be important biomarkers for breast cancer prognosis. Further studies with a larger sample size are warrant to further elucidate the association between TILs and breast cancer outcomes. in and can also be used in other tissue assays. Results from our pilot study indicated that increased iCD56 + cells were associated with worse OS of breast cancer patients. Increased sPD-L1 + cells and high sCD8+/FOXP3 + ratio were associated with more favorable OS in TNM stage III-IV breast cancer patients, whereas increased sFOXP3 + and total FOXP3 + cells indicated better OS among grade III breast cancer patients. CD8+/FOXP3 + ratio may be a better prognostic biomarker than CD8 or FOXP3 alone.


Introduction
Tumor-in ltrating lymphocytes (TILs), an important component of the tumor microenvironment (TME), represent a local immune response against cancer. Recent studies indicate that TILs have predictive value for immunotherapy and prognostic roles for various types of cancer [1]. Breast cancer was initially considered a poorly immunogenic tumor type; however, recent research has shown that some breast cancers have high levels of TILs, which have signi cant prognostic and predictive value, especially in human epidermal growth factor receptor 2 (HER2)-positive and triple negative breast cancers (TNBC) [2][3][4]. In solid tumors, TILs include T lymphocytes, macrophages, B cells, natural killer (NK) cells, and other immune cell subsets dispersed into both stromal and intratumoral compartments with varying spatial distributions [5]. Stromal TILs (sTILs) are in the brous stroma adjacent to tumor cells, while intratumoral TILs (iTILs) have cell-to-cell contact with tumor cells. The spatial interactions between TILs and cancer cells generate complex ecological dynamics that can ultimately impact tumor progression and response to treatment. Therefore, analysis of the spatial distribution and density of different immune cell populations that identify immune contexture components may be bene cial to cancer patients [6]. A recent study reported that high immune spatial scores, rather than immune cell abundance scores, were associated with poor recurrence-free survival in estrogen receptor (ER) + breast cancer [7]. The association between high spatial scores and late recurrence suggests a lasting memory of protumor immunity that may impact disease progression.
Hematoxylin and eosin (H&E)-based TIL assessment is easy for clinical application and has been proposed as a biomarker for inclusion in routine histopathological reporting to predict the prognosis or treatment response for breast cancer patients [8,9]. However, H&E-based TIL assessment is limited because it has a large inter-observer disagreement [10,11], cannot evaluate iTILs which are less frequent and more di cult to score than sTILs, and cannot differentiate various TIL subpopulations. Conventional immunohistochemistry (IHC) is a commonly used histopathologic technique that can examine TIL subsets and immune cell markers, such as cytotoxic CD8 + T cells, T regulatory cells (Tregs), and PD-L1 [12][13][14]. However, it has certain critical limitations, such as only labeling one marker per tissue section, high inter-observer variability in assessment, and lack of high-throughput quantitative analysis [15]. In recent years, various multiplex IHC (mIHC) techniques, such as chromogenic mIHC [16], cyclic immuno uorescent IHC [17], tyramide signal ampli cation (TSA)-based mIHC [18][19][20], metal-based mIHC [21,22], and DNA barcoding-based mIHC [23][24][25][26], have emerged to circumvent the limitations of conventional IHC. These techniques allow for the simultaneous detection of multiple markers on a single tissue section to provide comprehensive cellular spatial information and greater insight into the pathogenesis of cancer and responsiveness to immunotherapy. However, each technique has certain constraints [15], such as being time-consuming, costly, having low sensitivity, and a lack of standardized analysis software, making it impractical for most lab settings. Thus, it is imperative to develop an automated high-throughput methodology that is applicable in regular laboratory settings, feasible for large population studies, and scalable for the elucidation of clinicopathological signi cance of spatial TIL subpopulations on cancer outcomes.
In this study we established reliable and sensitive 5-plex uorescent mIHC methods in our regular IHC settings and developed high throughput automated quantitative analysis techniques for the automated evaluation of the density and spatial distribution of different TIL subsets. We then applied our approach in a pilot study to test associations between TILs and clinicopathological parameters and breast cancer prognosis.

Study participants
The human breast cancer samples evaluated in this study were from Nashville Breast Health Study (NBHS), a case-control study conducted in Nashville, Tennessee [27,28]. The NBHS included 2,694 breast cancer patients diagnosed between February 1, 2001 and December 31, 2011 that were identi ed through the Tennessee State Cancer Registry and ve major hospitals providing medical care for breast cancer patients. Participants were between the ages of 25 and 75, had no prior history of cancer other than nonmelanoma skin cancer, had a residential telephone, and spoke English. For breast cancer cases, information on ER, progesterone receptor (PR), and HER2 status of breast tumors was obtained from medical records or measured by the Vanderbilt Molecular Epidemiology Core Lab. Approval for this study was obtained from the institutional review boards of Vanderbilt University Medical Center and the other participating institutions, and all participants provided informed consent prior to enrollment in this study.
Tissue microarray (TMA) From a total of 1,255 tumor tissue samples collected by the NBHS, we randomly selected four TMA blocks, which include 208 patient samples for evaluation in this pilot study. Before TMA construction, all tissue samples were reviewed by an experienced study pathologist to con rm diagnosis. From each donor block, representative cancer areas of one spot in a peripheral region and two spots in central and middle regions were punched with a 1 mm needle and transferred into a recipient TMA block using a semi-automated Tissue Microarrayer. Each TMA block contained 170 tissue cores (13X13), including 52 cases of breast cancer, 13 internal control tissues/cell lines, and one orientation tissue core. The TMA slides were sectioned at 5 μm, covered with a thin layer of para n, and stored in a vacuum cabinet at 4 ºC until use.

Fluorescent multiplex immunohistochemistry (mIHC)
We developed two protocols for 5-plex uorescent mIHC to evaluate both TIL subpopulations and their spatial localization in breast tumor tissues. The rst mIHC panel included CD3 (for total T cells), CD8 (for cytotoxic T cells), and forkhead box protein 3 (FOXP3, for Tregs). The second panel included CD20 (for B cells), CD56 (for NK cells), and the immune check-point molecule programmed death ligand 1 (PD-L1).
Additionally, DAPI (for nuclei counterstain) and pan-cytokeratin (for tissue segmentation) were included in both panels. The mIHC methods were established based on PerkinElmer Opal TM Multiplex technology [18,29], optimized and validated by comparing conventional standard IHC (DAKO EnVision IHC kits) and single uorescence staining of each biomarker to produce speci c, non-interfering, and balanced labeling for three immune cell markers of interest per slide. All primary antibodies are commercially available and well validated [19,20,30,31]. A control TMA including 12 cores of human cancer and normal tissues was used for optimizing mIHC staining protocols and was stained in parallel with study TMA samples for quality control. The detailed mIHC staining methods are described in supplementary materials. After antigen retrieval followed by non-speci c blockings, slides were subsequently incubated with the rst primary antibody for 1 hour, Polymer-HRP Goat Anti-Rabbit/Mouse secondary antibody for 30 minutes, followed by TSA-Cy3 for ve minutes. After the microwave treatment for stripping antibodies followed by blocking steps, the slides were incubated with the second primary antibody at 4 ºC overnight. Then, slides were incubated with the secondary antibody as above, followed by TSA-FITC for ve minutes. The microwave treatment and blocking steps were repeated as above and the slides subsequently incubated with the third primary antibody for one hour, the secondary antibody for 30 minutes, and TSA-Cy5 for ve minutes. After microwave treatment and blocking steps, the slides were incubated with mouse anti-pan cytokeratin at 4 ºC overnight, followed by Cy7® goat anti-mouse IgG for one hour. Slides were thoroughly washed in PBS between each step. The stained slides were mounted with ProLong Gold Antifade Mountant with DAPI and stored in a box at 4ºC.

Automated image quanti cation
The mIHC stained slides were imaged using an Olympus IX-81 uorescence fully motorized and automated multispectral slide analysis system. Three cores per case were imaged with a low power objective (10X). Before applying the rst primary antibody, a pre-scan was performed to capture nonspeci c auto-uorescence signals in the FFPE tissues, which were then subtracted from speci c signals of stained slides. Grayscale images from the 5-plex mIHC were processed automatically using an established pipeline with the open-source imaging software CellPro ler v3.1.5 [32] to overlay each biomarker with segmentation markers (pan-cytokeratin and DAPI) for RGB images (16 bit tiff format). Based on proposed guidelines for semi-quanti cation of stromal TILs in breast cancer [8], TILs around ductal carcinoma in situ and normal lobules and those in tumor zones with crush artifacts, necrosis, and regressive hyalinization were excluded using the open-source ImageJ software. We then used CellPro ler to establish high-throughput image quanti cation pipelines for each biomarker. TILs within the tumor cell nest (intratumoral score for iTILs) and stromal area (stromal score for sTILs) were evaluated separately. TILs within 16 μm (45 pixels in our imaging system) to cancer cells were counted as iTILs. TIL intensity was recorded as positive count (the percentage of intratumoral or stromal TIL cell count over total cell count in tumor or stromal area) and positive density (the number of TILs per 100 μm 2 tumor or stroma cell nuclear area). Therefore, in addition to a total quantitative score for each biomarker, the spatial distribution of immune cell phenotypes in both intratumoral and stromal regions were also quanti ed.

Statistical analysis
Mean and standard deviation were applied to describe distribution for continuous variables and percent of each categories for categorized variables. Correlations between immune cell markers in human breast cancer tissue were examined using Spearman correlation coe cients. Cox proportional hazards models were employed to evaluate associations with overall survival (OS) for clinicopathologic characteristics of breast cancer and selected biomarkers. Tests for trends were conducted by entering the order of categorized variables as continuous variables in regression models. Multivariable Cox models included adjustment for tumor-node-metastasis (TNM) stage, tumor grade, ER and PR status, and race. All analyses were performed using Statistical Analysis Software (SAS O ce Analytics, version 14.3; SAS Institute Inc., Cary, NC) with a two-sided signi cance threshold of P < 0.05.

Establishment of uorescent mIHC and automated image quanti cation
We successfully stained and illustrated two panels of 5-plex immune cell markers (Figure 1). Staining protocols were optimized and validated by comparing conventional standard IHC (DAKO EnVision IHC kits) and single uorescent staining for each biomarker to produce speci c, non-interfering, and balanced labeling for three immune cell markers of interest per slide (Supplementary Figure 1). After imaging stained TMA slides, positive cells within the intra-tumoral area (iCD3+, iCD8+, iFOXP3+, iCD20+, iCD56+, and iPD-L1+), stromal area (sCD3+, sCD8+, sFOXP3+, sCD20+, sCD56+, and sPD-L1+), and total area of breast cancer tissues were scored separately by high-throughput automated quanti cation with CellPro ler. For each marker, TMA control tissues and sample cancer tissues were used as a training set for the optimization of measurement pipelines. For example, the automated imaging analysis of PD-L1 ( Figure 2) showed dependable results of spatial PD-L1 positivity in both intratumoral and stromal regions. Veri cation of all pipelines (Supplementary Figure 2) enabled a large batch of images to be automatically analyzed at one time. Assessment of biomarkers was conducted for three cores per case and then averaged for further statistical analysis.
Association of clinicopathological parameters with prognosis of breast cancer patients Among four TMAs, a total of 188 (89.4%) cases had evaluable tissue cores and were included in our analysis ( Table 1)

Association of TIL subsets and spatial expression with breast cancer prognosis
We evaluated the association of each individual biomarker as well as the CD8+/FOXP3+ ratio with overall survival (OS) of breast cancer patients ( Table 3). The results showed that higher iCD56+ cells and total CD56+ cells were associated with worse OS. Hazard ratios (HRs) for total mortality, comparing highest to lowest tertiles, were 4.89 (95% con dence interval (CI): 1.26-18.95) for iCD56+ cell density and 3.79 (95% CI: 1.12-12.87) for total CD56+ cell density. No other biomarkers had signi cant associations with OS when analyses included all breast cancer patients. TIL quanti cation using positive cell count percentage (Supplementary Table 1) and positive cell density (Table 3) showed similar results.
Associations with OS were further evaluated in analyses strati ed by race, TNM stage, grade, and ER/PR status (Table 4). There were no association between any of the six biomarkers and OS in subgroups strati ed by race (P>0.05). Both increased sPD-L1+ and sCD8+/FOXP3+ ratio showed more favorable outcomes in TNM stage III-IV subgroups. Breast cancer patients with the highest tertile of sPD-L1+ cells had a lower risk of mortality (HR, 0.01; 95% CI, 0.00-0.89; P=0.042 for trend) than those with the lowest tertile. Similarly, breast cancer patients with the highest tertile of sCD8+/FOXP3+ ratio had a lower risk of mortality compared with the lowest tertile (HR, 0.00; 95% CI, 0.00-0.12; P=0.036 for trend) in TNM stage III-IV patients. Neither CD8+ nor FOXP3+ cells showed signi cant association with OS in subgroups strati ed by TNM stage. Strati cation by tumor grade resulted in a trend of favorable prognosis with increased FOXP3+ cells among grade III breast cancer cases, with HRs for total mortality of 0.17 (95% CI: 0.03-0.97, P=0.037) for sFOXP3+ cell density and 0.07 (95% CI: 0.00-1.51, P=0.044) for total FOXP3+ cell density. TIL quanti cation using positive cell count percentage (Supplementary Table 2) and positive cell density (Table 4) showed similar results.

Discussion
The assessment of TIL subsets using mIHC and computer-assisted imaging quantitative analysis provides a robust tool for the accurate identi cation of TIL subsets and objective quanti cation of the immune context in tumors. In this study, the uorescent mIHC methods were established based on PerkinElmer Opal™ Multiplex technology, which was developed speci cally for FFPE tissue and showed high sensitivity in detection of low-abundance targets due to tyramide signal ampli cation [18,29]. Our three-day uorescent mIHC staining protocols can be easily conducted, even in resource-limited settings. After stained TMA sections were imaged, we used free and open source image analysis software for image editing and high-throughput analysis. First, the grayscale images from 5-plex mIHC were processed automatically by a pipeline established with CellPro ler [32] to overlay each biomarker with segmentation markers and to export as RGB images. Second, based on international TILs assessment guidelines [8], the unwanted areas were excluded using ImageJ. Third, for high-throughput automatic imaging quantitative analysis, we developed measurement pipelines for each biomarker using CellPro ler. The pipeline automatically quanti ed TILs in both the intratumoral (iTILs score) and stromal areas (sTILs score). The potential nonspeci c staining signals could be excluded by relating biomarker signal with nuclei. The individual TIL subset was quanti ed using both positive cell count and positive cell density. We found that both measures showed similar associations with breast cancer survival in this study. Hundreds of images were processed at one time, and quantitative data were exported automatically for statistical analysis. Hence, we successfully established a feasible computer-assisted high throughput automatic analysis method to evaluate the spatial distribution of TIL subpopulations in breast cancer, which is useful for further large population-based studies and for other in situ biomarker studies.
Although TILs have been widely observed in breast tumors, the prognostic potential of these markers and their predictive roles have only been investigated in the last decade [5][6][7]. Increasing evidence suggests a linear relation between high TIL levels and improved outcome in the patients with TNBC and HER2positive breast cancer [1,33,34] but an inverse association with survival outcomes in patients with ER + breast cancer [7,33]. TILs were also considered biomarkers for pathological complete response (pCR) in breast cancer patients treated with neoadjuvant therapy, with the positive correlation between pCR rate and TIL level [35,36]. The prognostic value of TILs, especially cytotoxic CD8 + T cells or immunosuppressive FOXP3 + Tregs, have been previously investigated in breast cancer; however, results remain inconsistent across different subtypes of breast cancer. In the current study, we found that the spatial distribution of TIL subsets was associated with overall survival among breast cancer patients (i.e., increased iCD56 + cells within the tumor bed area was associated with worse OS). Strati ed analysis indicated that higher expression of sPD-L1 and higher sCD8+/sFOXP3 + ratio correlated more favorably with OS in TNM stage III-IV breast cancer patients, whereas increased sFOXP3 + and total FOXP3 + cells were associated with better OS among grade III breast cancer patients. Furthermore, three T cell markers (CD3, CD8, and FOXP3) signi cantly correlated with each other, and PD-L1 positivity signi cantly correlated with all other markers, especially CD3+, CD8+, and CD56 + cells in stromal and intratumoral areas ( Table 2). These ndings are generally consistent with the potential biological impact of selected biomarkers on breast cancer cells, as described below.
CD56 + NK cells are historically considered as the rst line of host defense against tumor cells, and in ltration of tumors with NK cells is a prognostic marker for several malignancies, including breast cancer. A recent study showed that the frequency of NK cells increased signi cantly in poorly differentiated (grade III) breast cancer and in tumor draining lymph nodes, suggesting a suppressive phenotype for these cells in breast cancer [37]. Another study showed CD56 + NK cell in ltration was inversely correlated with PR and ER receptor expression status (p = 0.021 and p = 0.007, respectively), two well-established biomarkers of better prognosis in breast cancer [38]. In agreement with those ndings, our study found that total CD56 + NK cells and iCD56 + cell density was inversely associated with OS.
There are two major NK subpopulations that have been described in humans: a CD56 dim subset that predominates in blood and exhibits a high cytotoxic potential, and a CD56 bright subset that is more abundant in secondary lymphoid tissues and functions as regulatory and tolerant subsets [39]. The distinct suppressor subsets of NK cells in the TME could suppress antitumor immune responses in solid tumors through inhibitory receptors on the surface of NK cells to downregulate activating signals, leading to a decline of the antitumor immune response and resulting in poor outcomes for cancer patients. The CD56 + NK cells detected in this study might belong to the subset of immunosuppressive NK cell phenotypes.
Previous pooled studies indicate that CD8 + T cells are associated with better disease-free survival (DFS) and breast cancer speci c survival (BCSS) [2]. A recent meta-analysis with 37 studies also found that a higher CD8 + TIL level was associated with better DFS, although no signi cant association was found with OS among TNBC patients [40], which is largely consistent with the results of our study. The prognostic role of Tregs, de ned as FOXP3 + T cells, remains controversial [12,13]. A meta-analysis of twenty-ve published studies with 22,964 patients reported that high levels of Treg lymphocytes were associated with poor DFS and OS, but not with poor BCSS [2]. In contrast, Bottai et al. reported that FOXP3 + cells were signi cantly associated with better RFS and OS among early stage TNBC patients. However, the prognostic value of FOXP3 + cells was not signi cant after adjusting for CD8+ [41]. An updated meta-analysis of thirty-seven studies with 10,258 patients showed that higher FOXP3 + TIL level was associated with better DFS but not OS in TNBC patients [40]. Furthermore, a higher CD8/FOXP3 ratio was found to be related to higher pCR rate and improved DFS and OS among advanced stage breast cancer patients [42]. In our study, we saw positive associations between stromal FOXP3 + cells and the CD8+/FOXP3 + ratio with better OS, especially in grade III or stage III-VI patients. Our results suggest that the CD8+/FOXP3 + ratio may be a better prognostic biomarker than CD8 or FOXP3 alone for advanced stage or high-grade breast cancer.
Immune check-point molecules programmed cell death 1 (PD-1) and PD-L1 are key immune response modi ers through regulating T-cell activation and immune surveillance [43]. Expression of PD-L1 is not only related to the response of immune checkpoint therapy but also correlated with prognosis in many cancer types [44]. A recent meta-analysis with a total of 14,367 primary breast cancer patients showed that PD-L1 expression on tumor cells was associated with shorter DFS and OS, while PD-L1 + TILs was related to better DFS and OS [45]. In this study, we found that higher sPD-L1 + cells correlated more favorably with OS among TNM stage III-IV breast cancer patients, indicating that PD-L1 + TILs may be an indicator for favorable prognosis in breast cancer patients with advanced stage disease.
Our study has several strengths. Our newly established reliable and sensitive uorescent mIHC and automatic computer-assisted quanti cation methods using free and open source imaging software are easy to implement, even in resource-limited settings. The mIHC work ow provided a cost-e cient complete solution from staining to imaging to analysis protocol and can facilitate high throughput analysis that is suitable for large clinical trials or epidemiologic studies. The positive cell percentage and density of TIL markers in intratumoral and stromal compartments were quanti ed automatically and consistently through whole study samples. In addition, the results generated from analysis of our TMA sections are consistent with the potential biological roles of our selected biomarkers on breast cancer cells and on the prognosis of breast cancer patients, which also supports that three 1 mm cores for each breast cancer case for TIL and other biomarker studies is su cient [46][47][48][49][50][51][52][53].
Our study also has some limitations. Information on cancer recurrence and cause of death was not available in the NBHS cohort study, so we were unable to conduct disease-free and breast cancer-speci c survival analysis. The small sample size of this pilot study meant low statistical power to detect moderate associations between TILs and breast cancer outcomes and prohibited additional analyses by breast cancer subtype. Further studies with larger populations are needed to elucidate the association between TILs and breast cancer outcomes.
In conclusion, we developed a high throughput work ow for uorescent mIHC and automatic quanti cation analysis that enables investigation of the expression and spatial distribution of different TIL subsets and immune biomarkers in large scale studies and can also be used in other tissue biomarker assays. Results from our pilot study indicated that increased iCD56 + cells were associated with worse OS of breast cancer patients. Increased sPD-L1 + cells and high sCD8+/FOXP3 + ratio were associated with more favorable OS in TNM stage III-IV breast cancer patients, whereas increased sFOXP3 + and total FOXP3 + cells indicated better OS among grade III breast cancer patients. CD8+/FOXP3 + ratio may be a better prognostic biomarker than CD8 or FOXP3 alone. Abbreviations BCSS, breast cancer speci c survival; CI, con dence interval; DFS, disease-free survival; ER, estrogen receptor; FOXP3, forkhead box protein 3; H&E, hematoxylin and eosin; HER2, human epidermal growth factor receptor 2; HR, hazard ratios; mIHC, multiplex immunohistochemistry; NBHS, Nashville Breast Health Study; NK, natural killer cells; OS, overall survival; pCR, pathological complete response; PD-L1, programmed death ligand 1; PR, progesterone receptor; TILs, tumor in ltrating lymphocytes; iTILs, intratumoral TILs; sTIL, stromal TILs; TMA, tissue microarray; TME, tumor microenvironment; TNBC, triple negative breast cancer; TNM, tumor-node-metastasis; Tregs, T regulatory cells; TSA, tyramide signal ampli cation.

Declarations
Acknowledgements Dr. Robert Coffey's lab at VUMC kindly provided technical support for the establishment of 5-plex urescent mIHC and TMA image capturing. We also thank Ms. Nabela Hamm for tissue processing of the NBHS samples, Dr. Mary Shannon Byers and Ms. Rachel Mullen for assistance in editing and submitting this manuscript. Funding NBHS is supported by the NIH grant R01CA100374. NHBS Data and tissue sample collection, TMA construction, uorescent mIHC staining, and imaging quantitative analysis were conducted in the Survey and Biospecimen Shared Resource, which is supported in part by the Vanderbilt-Ingram Cancer Center (P30CA68485).
Data availability The authors con rm the data that have been used in this work are available on reasonable request.
Compliance with ethical standards