Automated Analysis of Lymphocytic Infiltration, Tumor Budding, and Their Spatial Relationship Improves Prognostic Accuracy in Colorectal Cancer

Tumor–node–metastasis (TNM) staging remains the gold standard for stratification of colorectal cancer patients into prognostic subgroups and is fundamental for treatment selection (1, 2). Surgical resection without the use of adjuvant treatment results in long-term disease-free survival in most stage II colorectal cancer cases. However, 20% of patients experience disease-specific death (3). In addition, some stage III patients experience better outcomes than some stage II patients. This is in spite of revisions to the TNM guidelines, which, for example, subdivide pT stage II into three categories (4). Alongside pT stage, tumor differentiation aids the pathologist in substratifying patients into high or low risk of tumor progression; however, grade is highly subjective (5). Tumor budding, which morphologically resembles small, poorly differentiated clusters of cells (1–4 cells) disseminating at the invasive margin (IM), has consistently been correlated with poor colorectal cancer prognosis (6–8).

Both TNM staging and tumor budding concentrate solely on the tumor, but it is now well established that the host interaction, such as immune infiltrate, within the tumor microenvironment (TME) plays a crucial role in tumor progression (9, 10). Specifically, in colorectal cancer, the pioneering studies of Galon and colleagues demonstrated that CD3⁺ and CD8⁺ T-cell infiltration was an independent prognostic factor predicting patient survival (9, 11). Tumor buds (TBs) are hypothesized to represent a highly invasive tumor subpopulation, especially if they also have the ability to evade the host immune response (6–8). The interaction between immune response and TBs may therefore hold vital information pertaining to the tumor's potential to metastasize.

Previous research has demonstrated the ability to automatically quantify CD3⁺ and CD8⁺ T cells within the TME, thereby achieving a greater level of standardization than with subjective, manual reporting (9). Similarly, we have previously shown that tumor budding can be standardized through automated image analysis (12, 13). Although the automatic quantification of both immune infiltrate and tumor budding has been associated with patient prognosis, they have not been studied together. Furthermore, by coregistration on the same tissue section, it is possible to accurately study the spatial interaction with each other. We apply such an approach in this study to investigate if quantifying both features alongside their spatial interaction allows for a more precise and accurate prognosis for stage II colorectal cancer patients.

This study included three independent cohorts of stage II colorectal cancer patients who underwent surgical resection in Edinburgh, UK, hospitals over the years 2002–2004 (n = 170) and the National Defense Medical College Hospital (NDMCH), Japan, over the years 2006–2011 (n = 62). The formalin-fixed paraffin-embedded (FFPE) block that contained the deepest cancer invasion for each patient was selected after the review of hematoxylin and eosin (H&E) slides. Cancer blocks from Edinburgh hospitals comprised specimens of cross-sectional cut, whereas the ones from the NDMCH were of longitudinal cut. Associated clinical data included features from the original pathology report, such as TNM staging and differentiation as well as detailed follow-up information including disease-specific survival (DSS): up to 11.5 years for the Edinburgh cohort and 8.6 years for the Japanese cohort. One hundred fourteen consecutive patients who underwent surgical resection in Edinburgh hospitals over the years 2002–2003 were used as a training cohort for this study. All available patients who underwent surgical resection in Edinburgh hospitals in the following year (2004) were used as a first validation cohort (n = 56). We then used a similar number of patients who underwent surgical resection in the NDMCH, Japan, over the years 2006–2011 as an international validation set (n = 62). The clinicopathologic data of all the cohorts are shown in Table 1. This study was conducted in accordance with the Declaration of Helsinki. Ethical approval was obtained after review by the NHS Lothian NRS BioResource, REC-approved Research Tissue Bank (REC approval ref: 13/ES/0126), granted by East of Scotland Research Ethics Service and by the Ethics Committee of the National Defense Medical College (approval ref: No. 2992). The ethical approval permitted us to use archival diagnostic samples once they were deidentified for the purposes of the research.

We selected externally validated specific primary antibodies: primary antibody against CD8⁺ cytotoxic T cells (monoclonal mouse anti-human CD8, M7103; Dako) was previously used for the international validation of the consensus IMMUNOSCORE®* by Pages and colleagues (14), primary antibody against pan T cells (polyclonal rabbit anti-human CD3, A045201-2, Dako) used by Harter and colleagues (15), in assessment of tumor-infiltrating lymphocytes in human brain metastases and for pan-cytokeratin (PanCK; mouse monoclonal anti-human cytokeratin, M351501-2, Dako) previously used by Koelzer and colleagues (16), for the assessment of TBs in colorectal cancer.

In order to assess the in-house antibody performance, each antibody was individually assessed by brightfield IHC. The primary antibodies against CD3 and CD8 T cells were each tested individually on 3-μm thick serial sections of normal tonsil FFPE tissue, known to have high lymphocytic cell concentration. For assessment of PanCK primary antibody, a colorectal cancer optimizing tissue microarray, consisting of 100 cores taken from primary tissue blocks, was used. Expression of all cell markers was detected using 3,3′-diaminobenzidine (DAB) chromogen (K3468; Dako). Nuclei were counterstained with hematoxylin (RBA-4213-00A; CellPath). In parallel, H&E (RBB-0100-00A; CellPath) staining was performed on serial sections and specificity of antibodies from the IHC was assessed and compared visually by an experienced pathologist (D.J. Harrison) and research scientists (I.P. Nearchou and P.D. Caie). No CD8-positive but CD3-negative cells were present, giving confidence to the antibodies' specificity. Positive and negative controls were used in each IHC staining; human tonsil FFPE tissue without primary antibody for CD3 and CD8 T cells and a colorectal cancer optimizing tissue microarray without primary antibody for PanCK. Western blot analysis was also used to assess the specificity of the PanCK antibody using colorectal cancer cell lines.

Once specificity of antibodies was confirmed in brightfield, each target was assessed by a uniplex immunofluorescence. Uniplex immunofluorescence was performed on serial sections of the colorectal cancer optimizing tissue microarray (specified above) with the use of individual tyramide signal amplification (TSA) fluorescence kits. After deparaffinization, slides were placed in a microwaveable pressure cooker containing boiling antigen retrieval (AR) buffer (0.1 M sodium citrate pH6) and microwaved for 5 minutes at 75°C. Slides were allowed to cool down in running tap water for 20 minutes at room temperature and were then rinsed with 0.1% phosphate-buffered saline-Tween 20 (PBS-T). PBS-T was used for all washing steps. Tissues were then blocked with 3% hydrogen peroxide for 5 minutes and Dako serum-free protein block (X090930-2, Dako) was then applied for 10 minutes. Slides were then incubated for 30 minutes with the same antibodies as for IHC at these dilutions: CD3 (dilution 1:400), CD8 (dilution 1:200), and PanCK (dilution 1:100). Slides were then washed (3 × 5 minutes). The slides that had previously been incubated with CD3 or CD8 antibodies were then incubated in prediluted anti-mouse or anti-rabbit HRP-conjugated secondary antibodies (K400111-2 and K400311-2, respectively; Dako) for 30 minutes. Following three more washes of 5 minutes, slides were incubated for 10 minutes with either TSA fluorescein (FITC) for CD3 antibody visualization or TSA Cyanine 5 for CD8 antibody visualization (NEL741B001KT and NEL745B001KT, respectively; PerkinElmer) at 1:100 dilution. The slide that had been previously incubated with PanCK antibody was then incubated in anti-mouse AlexaFluor 555 (A21422, Thermo Fisher Scientific) in antibody diluent (S080983-2, Dako) at 1:100 dilution for 30 minutes. Finally, following a washing step, the slides were counterstained for 10 minutes with Hoechst 33342 (H3570, Thermo Fisher Scientific) at 1:20 dilution in deionized water. Sections were then dehydrated in 80% ethanol and mounted with ProLong Gold Antifade Reagent (P36930, Thermo Fisher Scientific). Optimal antibody dilutions in the uniplex immunofluorescence slides were chosen to ensure best signal to noise and specific staining pattern when digitized on the Zeiss Axioscan.z1 whole slide scanner (Zeiss). Fluorescence slides were compared with serial sections of brightfield DAB slides using the same antibodies. Positive and negative controls were included in each staining run. The positive controls for the CD3 and CD8 antibodies were the tonsil tissue as specified above.

Once the appropriate dilution of primary antibodies, secondary antibodies, and fluorophores were found under uniplex immunofluorescence, we then assessed them multiplexed on a single tissue section. In order to eliminate species cross reactivity, microwave stripping was performed between the CD8 and PanCK staining. Microwave stripping was performed as follows. AR buffer was heated in a pressure cooker in the microwave for 12 minutes prior to adding the slides in the solution. Slides were then microwaved for 17 minutes using the autodefrost function. Sections were then allowed to cool down in running tap water for 20 minutes and washed for 5 minutes. The immunofluorescence protocol utilized to label the study samples was performed as above using a Dako Link48 Autostainer (Dako). The first run of labeling included the primary antibody against CD8 prior to antibody stripping, and labeling with a primary antibody cocktail of CD3 and PanCK prior to counterstaining with Hoechst. The multiplexed section was compared with the uniplex labeled serial sections and we confirmed that there was no increase or decrease in the level of specific signals or nonspecific background signals.

Whole slide fluorescence images were captured with a 20× objective using an Axioscan.Z1 (Zeiss). Exposure times for Hoechst, Cy3, Cy5, and FITC channels were 12, 120, 8, and 25 ms, respectively.

Digital whole slide fluorescence images, in Zeiss' .czi file format, were uploaded into HALO Next-Generation Image Analysis software (version 2.1.1637.18; IndicaLabs) for image analysis. Full algorithm workflow and settings are shown below:

All nuclei in the whole slide image were automatically segmented using the commercially available High-Plex FL (version 2.0) module within the HALO Next-Generation Image Analysis software. Optimized module parameters for nuclei detection included dye weights (Hoechst nucleus weight = 4, Cy5 nucleus weight = 0.232, FITC nucleus weight = 0.232), nuclear contrast threshold (0.5), minimum nuclear intensity (0.079), nuclear segmentation aggressiveness (0.65), and using the default nuclear size setting (1–500 μm²). All samples were analyzed using these specific settings.

Within the same software module in which we segmented the nuclei (High-Plex FL), cells were then classified as CD3 (FITC) or CD8 (Cy5) positive based on dye nucleus-positive threshold (Cy5 = 0.132, FITC = 0.15), cytoplasm-positive threshold (Cy5 = 0.075, FITC = 0.5) and membrane-positive threshold (Cy5 = 0.075, FITC = 0.5). The algorithm settings were kept constant across all patients. Next, the algorithm automatically quantified the number of each lymphocyte classification across the whole slide image. The flood tool was used to identify the invasive front. An IM was created that consisted of a width span of 500 μm in and out from the invasive front. The tumor core area (CT) was set to be the whole tumor area excluding the 500 μm before the invasive front (Fig. 1A). Lymphocytes (CD3⁺ and CD8⁺ T cells) infiltrating both the CT and the IM of the tumor were automatically quantified (Fig. 1B) and their densities (cells/mm²) were exported to assess their prognostic significance.

LIS was performed based on the methodology described previously by Galon and colleagues (9). Using the lymphocytic infiltration analysis results, CD3⁺ and CD8⁺ lymphocytes were counted in both the IM and CT, resulting in the reporting of the following features: CD3⁺ T cells in the CT, CD3⁺ T cells in the IM, CD8⁺ cells in the CT, and CD8⁺ T cells in the IM. Optimal cutoff points for the four features were calculated based on survival data of the training cohort (Supplementary Table S1). Values above the optimal cutoff points for each feature were scored 1, whereas values below were scored 0. A LIS was then calculated by the sum of all scores for each patient sample. Finally, similarly to Liu and colleagues (17), patients were stratified into two groups: high LIS > 2 and low LIS ≤ 2.

The PanCK intensity within each whole slide image was initially assessed using the Area Quantification FL v1.2 module without any changes to the commercially available module thresholds. Patient samples were grouped into two categories based on their resulting PanCK intensity: low (average dye intensity < 2.16 × 10⁻²) and high (average dye intensity > 2.16 × 10⁻²). Using samples from within each PanCK group (high or low), machine learning classifiers were trained to segment tumor from nontumor regions. Two High-Plex FL (version 2.0) modules were next created: one using the classifier trained from low PanCK patient samples and one from the high PanCK patient samples. Patient samples with an average PanCK intensity below 2.16 × 10⁻² were analyzed to segment tumor from nontumor using the classifier termed “low PanCK” and the other samples with the classifier termed “high PanCK.” These modules were run within a 1,000-μm border inward from the invasive front, termed the tumor-budding region of interest (TBROI; Fig. 1C). Within the tumor regions detected, nuclei were detected using the methodology described above, and cancer cells were classified based on PanCK positivity in nucleus (0.1), cytoplasm (0.08), and membrane (0.05; Fig. 1D). These cancer cell classification thresholds were kept constant across the entire cohort of images. TBs were classified as tumor clusters containing up to 4 PanCK⁺ cells. TB density (TBs/mm²) and number were then exported from the algorithm. Mucin pools and areas of necrosis were excluded.

Spatial coordinates of the lymphocytes and TBs were imported into a spatial plot within the HALO software. From this plot, the Spatial Analysis algorithm was utilized to calculate the number of CD3⁺ and CD8⁺ T cells within a range between 0 and 100 μm radii of TBs (Fig. 1E) in consecutively increasing 10 μm steps, therefore creating 11 classes: 0–10 μm, 0–20 μm, 0–30 μm, and so on. All classes of radius size were significantly associated with survival; however, after assessing the prognostic significance of iterative groupings of varying radius sizes, there was a statistical difference observed between the 0–50 μm and the 0–100 μm groupings. The classes were therefore merged into two categories: 0–50 μm and 0–100 μm radii of TBs. Finally, the number and densities of CD3⁺ and CD8⁺ T cells within these distances were established.

All features exported from the image analysis are listed in Supplementary Table S2.

All statistical analyses, unless otherwise stated, were performed in RStudio 1.1.419 running R 3.4.3 (18, 19). Continuous image analysis data, categorical clinical data, and survival data were loaded into R. The relationship between lymphocyte distribution patterns, TBs, and patient characteristics was assessed using the χ² test. For comparison with categorical patient data, image analysis features were divided into low and high groups according to optimal cutoff points based on survival data. To assess the relationship between lymphocytes and TBs, Spearman correlation coefficients (r) were applied on continuous data. All r = 0 values indicated no association between features, r > 0 indicated a positive association where as one variable increases, so does the other, and r < 0 indicated a negative association where as one variable increases the other decreases. Associations between features were assessed on the training cohort.

Univariate Cox regression was performed on both the categorical clinical and continuous image analysis data. All P < 0.05 values were considered statistically significant.

The least absolute shrinkage and selection operator (LASSO) penalized Cox proportional hazard regression was performed for identification of significant features (20), using the glmnet package (21). The significant features identified by LASSO were then inputted into a random forest (RF; n = 500) decision-tree model using the randomForest package (22) and were ranked according to their mean decrease Gini coefficient. Features with a mean decrease Gini of above 4 were taken forward for further analysis. Optimal cutoff points for these four features were calculated, and iterative combinations of these features were tested. The model with most predictive value was then selected, termed the “tumor bud-immuno spatial index” (TBISI). In order to avoid overfitting during the development of the model, validation was performed for both LASSO Cox regression model (10-fold cross-validation) and RF (out of the bag). Univariate Cox regression was performed for the TBISI, pT stage, TB number, and LIS, and a bootstrap resampling technique was applied (number of samples = 1,000) in SPSS 24 (23). Multivariate Cox regression with a forward stepwise method and the receiver operating characteristic (ROC) curve were also applied on to compare the TBISI, pT stage, TB number, and LIS in SPSS 24. The prognostic value of TBISI was then assessed on two independent validation cohorts. All validations were performed using the fixed optimal cutoff points of the training cohort. All optimal cutoff points for the categorization of image analysis features including tumor budding, the LIS, and for the development of the TBISI were based on DSS and were calculated using the survminer package (24). All Kaplan–Meier (KM) curves were plotted using the survival package (25). The P values from KM curves were corrected for false discovery rate (FDR) using the Benjamini–Hochberg procedure to account for multiple hypothesis testing (26).

Study workflow is illustrated in Supplementary Fig. S1.

The clinicopathologic characteristics of all three cohorts of stage II colorectal cancer patients are shown in Table 1. This study cohort comprised a training cohort of 114 stage II colorectal cancer patients, of which 57 were female and 57 were male. The validation cohorts included 34 male and 22 female patients in validation cohort 1 and 42 male and 20 female patients in validation cohort 2. The range of patients' age was 37–96 years old. There were 43 right-sided colon carcinomas, 38 left-sided, and 33 rectal cancers in the training cohort and 29 and 14 right-sided colon carcinomas, 9 and 22 left-sided, and 18 and 26 rectal cancers in validation cohorts 1 and 2, respectively.

Automated image analysis was performed to quantify TBs (within TBROI), CD3⁺ cell, and CD8⁺ T-cell densities (within IM and CT). Optimal cutoff points based on survival data were used to categorize image analysis features of the training cohort into low and high groups (Supplementary Table S1). TB number and density were highly correlated with each other (correlation coefficient = 0.7466; P < 2.2 × 10⁻¹⁶); therefore, TB number was selected to be used for correlation assessments. The relationship between patient characteristics and both lymphocyte distribution patterns and TB number was then assessed. A high density of CD3⁺ T cells at the IM was associated with no extramural lymphovascular invasion (EMLVI; P = 0.049), whereas in the CT, it was correlated with moderate differentiation status and nonmucinous histologic type (P = 0.007 and P = 0.023, respectively). A high density of CD8⁺ cells at the IM and across the whole tumor section (WTS) was significantly associated with female gender (P = 0.016 and P = 0.030, respectively). Finally, advanced pT stage was significantly associated with a high TB number at the TBROI (P = 0.031; Table 2).

TB number was tested for association with lymphocyte densities of the training cohort, using the Spearman correlation coefficient (Fig. 2). This indicated a weak negative relationship between TB number and all lymphocytic cell densities (r < 0) regardless of distribution pattern (i.e., IM, CT, WTS) or lymphocytic cell subpopulations (CD3⁺ or CD8⁺). However, a higher density of CD3⁺ T cells at the IM, CD8⁺ T cells at the IM, and WTS were significantly correlated with a lower TB number (P = 0.016, P = 0.037, and P = 0.041, respectively; Supplementary Table S3).

To assess the prognostic significance of the categorical clinicopathologic data of the training cohort, univariate Cox regression was applied. KM survival analysis was used to further assess any significant prognostic features from the clinicopathologic data of the training cohort (Fig. 3). Univariate Cox regression showed that only pT stage (HR = 4.143; 95% CI, 1.480–11.570; P = 0.006) was significantly associated with DSS (Table 1). KM survival function was then applied and showed that patients of T4 stage have poor DSS (73% survived) compared with T3 stage patients (91% survived; P = 0.003; Fig. 3A).

To explore the prognostic effects of the quantified lymphocytic distribution patterns, TB number and density within the training cohort, we performed univariate Cox regression analysis. This analysis was also performed for the spatial relationships between lymphocytes and TBs. The following specific continuous features were calculated and exported from the image analysis algorithm: (i) in IM, CT, and WTS; densities of CD3⁺ and CD8⁺ T cells and ratio of CD3⁺ to CD8⁺ T cells; (ii) TB number and density at the TBROI, (iii) within 0–50 μm and 0–100 μm radius of TBs; number and mean number of CD3⁺ and CD8⁺ T cells. The densities of both CD3⁺ and CD8⁺ T cells were significantly associated with DSS (P < 0.05). TB number and density conferred prognostic significance (HR = 1.001; 95% CI, 1.000–1.001; P = 0.0004 and HR = 1.030; 95% CI, 1.000–1.060; P = 0.0446, respectively). KM analysis showed TB number to be a significant prognostic factor; 95% of patients in the low tumor budding group survived compared with the high tumor budding group, in which 73% of patients survived (P = 0.0006; Fig. 3D). Furthermore, a high mean number of lymphocytes near TBs was significant in predicting DSS for stage II colorectal cancer patients (P < 0.05; Supplementary Table S2).

All clinicopathologic and image analysis data of the training cohort (listed in Supplementary Table S2) were input for a LASSO penalized Cox proportional hazard regression. This was performed to identify the features that add significant value to the prediction of DSS over time and therefore to reduce the data dimensionality. The analysis reported that there were eight significant candidate features (CD3⁺ T-cell density at WTS, mean CD3⁺CD8⁺ T-cell number within 0–50 μm of TBs, TB number, CD8⁺ T-cell density in CT, pT stage, EMLVI, age, and differentiation). The eight features that added value to the model were used as further input into a RF decision-tree model for optimal prognostic feature selection. The RF model reported four features with a mean decrease Gini of above 4 (CD3⁺ T-cell density in WTS, mean CD3⁺CD8⁺ T-cell number within 0–50 μm of TBs, TB number, CD8⁺ T-cell density in CT; Table 3).

Optimal cutoff points for these four features (CD3⁺ T-cell density in WTS, mean CD3⁺CD8⁺ T-cell number within 0–50 μm of TBs, TB number, CD8⁺ T-cell density in CT) were calculated for the training cohort using the patient survival data (389.6 cells/mm², 4.1, 1104.0, and 114.7 cells/mm², respectively). The least significant feature was then removed in an iterative process until the removal of a feature decreased the prognostic value of the model. This led to the creation of the TBISI consisting of three features: CD3⁺ T-cell density in WTS, mean CD3⁺CD8⁺ T-cell number within 0–50 μm of TBs, and TB number. Patients with CD3⁺ T-cell density in WTS and mean CD3⁺CD8⁺ T-cell number within 0–50 μm of TBs below the cutoff point and TB number above the cutoff point were classified as a “high risk” of disease-specific death group; the remainder of the patients were classified as a “low-risk” group. KM curves and bootstrap univariate Cox regression P values were used to determine the significance of this index. The created TBISI was highly significant in the prediction of disease-specific death from both KM analysis and Cox regression (HR = 18.75; 95% CI, 6.46–54.43; P = 1.202 × 10⁻¹³; Fig. 3J). KM and univariate Cox-regression analyses were also applied to pT stage (HR = 4.143; 95% CI, 1.48–11.57; P = 0.0033; Fig. 3A), TB number (HR = 5.899, 95% CI, 1.875–18.55; P = 0.0006; Fig. 3D) and LIS (HR = 7.02; 95% CI, 2.49–19.84; P = 1.922 × 10⁻⁰⁵; Fig. 3G). Multivariate Cox regression model (forward stepwise) revealed that this TBISI was the only feature that added value to the model when pT stage, TB number, LIS, and TBISI were included (Supplementary Table S4). Finally, the ROC curve showed that the TBISI had the greatest area under the curve (area = 0.785; 95% CI, 0.629–0.941) when compared with pT stage or each of the model's component individually (Supplementary Fig. S2).

The TBISI developed from the training set, alongside the LIS and TBs, were evaluated on two independent validation cohorts: validation cohort 1 patients from Edinburgh hospitals (n = 56) and international validation cohort 2 patients from the NDMCH, Japan (n = 62). Cutoff points calculated from the training set data were directly applied to the two validation cohorts. Univariate Cox regression and KM analyses were calculated using the patient data from the two validation cohorts and used to assess the clinicopathologic data, TB number, LIS, and the model TBISI (Table 1). pT stage was significantly associated with patient survival in validation cohort 1 (HR = 5.00; 95% CI, 1.52–16.46; P = 0.008; Fig. 3B), but not in the Japanese validation cohort 2 (Fig. 3C). TB number was significant in both validation cohorts 1 and 2 (HR = 5.29; 95% CI, 1.14–24.50; P = 0.033 and HR = 8.816; 95% CI, 2.70–28.82; P = 0.0003; Fig. 3E and F, respectively), whereas the LIS was only significant in validation cohort 1 (HR = 4.72; 95% CI, 1.02–21.87; P = 0.047) (Fig. 3H), although trended to significance in validation cohort 2 (Fig. 3I).

TBISI was highly significant on both validation cohorts; validation cohort 1 with HR = 12.27; 95% CI, 3.52–42.73; P = 8.2 × 10⁻⁰⁵ and in validation cohort 2 with HR = 15.61; 95% CI, 4.69–51.91; P = 7.42 × 10⁻⁰⁶ (Fig. 3K and L, respectively). The TBISI demonstrated an increase in hazard ratio and significance in both validation cohorts than when analyzing the sum of its parts independently.

Risk assessment and thus treatment choices for stage II colorectal cancer patients are in need of improvement. A growing awareness exists that the TME contributes to cancer differentiation, progression, invasion, and metastasis (27–29). In colorectal cancer, the two most promising candidate prognostic features involve two distinct aspects of the TME: tumor budding (2, 30–32) and the IMMUNOSCORE®* (9, 10, 14, 33). Tumor budding represents small clusters of infiltrating cancer cells within the IM, whereas the IMMUNOSCORE®* quantifies tumor-infiltrating lymphocytes (TILs) within both the IM and the CT. However, these features are predominantly assessed independently of each other. The question, therefore, arises as to whether tumor budding and lymphocytic infiltration combined augment their individual prognostic significance and whether their association and spatial relationship adds further prognostic value for stage II colorectal cancer patients. The work presented here utilizes automated image analysis to quantify and compare TBs and TILs coregistered to the same tissue section.

We report four major findings: first, we have shown significant associations between clinicopathologic factors, CD3⁺ CD8⁺ T-cell densities, and TB number; second, we have shown that LIS and tumor budding are strong independent prognostic factors for DSS in stage II colorectal cancer; third, we have demonstrated that the spatial interaction between lymphocytes and TBs holds additional prognostic significance; and finally, a standardized prognostic index was created combining all the above features to achieve a higher significant cohort stratification than any individual prognostic feature tested.

Emerging technologies in tissue analysis have allowed further understanding of the role of the TME in disease progression. Automated image analysis has been shown to have a high potential for unbiased and reproducible assessment of histopathologic features in colorectal cancer, such as for TBs, lymphatic vessel invasion (12), and TILs (14), as well as in other disease types (34, 35). To an extent, this is assisted by the advancement of acquiring multiparametric data from a single tissue section, which can provide a more complete understanding of the underlying mechanisms within the TME. Even though some multiplex IHC or immunofluorescence methodologies have shown the ability of high-level multiplexing on a single FFPE tissue sample, the analysis may be restricted to manually selected regions of interest or through their applications on tissue microarrays (36). Such methodologies can also be associated with long acquisition times for whole tissue sections and sequential stitching of image tiles into one large image (37). In this study, we quantified TILs and TBs coregistered to a single whole slide section through the use of four-plexed immunofluorescence, whole slide scanning, and automated image analysis. This methodology permitted the colocalization of the four tissue markers, required to quantify the features measured in this study, across the whole tissue section and at single-cell resolution. Spatial coordinates of each cell type were exported from the same tissue section ensuring more accurate analysis than is achieved from coregistering serial sections. This methodology therefore conserves precious tissue, reduces scanning time and cost, while also allowing the objective analysis of the entire TME and its cellular interactions.

Our results showed that the presence of high CD3⁺ T-cell density correlated with favorable clinicopathologic features, at the IM with the absence of EMLVI and at the CT with moderate differentiation status and nonmucinous histologic type. These data concur well with prognostic analyses on the propitious impact of overall TIL densities in colorectal cancer patients (9, 10). We also found that a high CD8 T-cell density was correlated with the female gender, a finding that has previously been shown in invasive colorectal cancer by Jang (38). However, the mechanisms underlying this observation remain unknown. An elevated number of TBs was significantly associated with advanced pT stage, correlating favorably with previous literature showing TBs being associated with adverse prognostic features (39–41). These findings highlight the impairing effect T cells have on tumor progression, while suggesting that high tumor budding cancers harbor a more aggressive phenotype.

Our data indicate a significant inverse relationship between the densities of TILs and TB number, which substantiates the work of Zlobec and colleagues (42), who found the absence of TILs to be associated with the presence of tumor budding. This may suggest that during the transition of TBs from epithelial-to-mesenchymal phenotype, they might occupy a niche that may not elicit an immune response. Indeed, it has been suggested that TBs lose MHC expression during this transition (7), and others have shown that TB positivity is significantly associated with PD-L1 positivity (43). The data presented here add further support to the concept of TBs having the ability to evade the antitumoral host response and thereby embark on the first steps of the metastatic cascade.

Despite repeated evidence of tumor budding conferring prognostic significance in colorectal cancer (8, 12, 41, 44), its routine clinical assessment has not yet been established due to the lack of a standardized methodology of assessment. Most studies investigating tumor budding rely upon the use of manual quantification, often based on subjective selection of high magnification regions of interest (2). The use of cytokeratin IHC has previously been shown to improve the detection of TBs as well as increasing interobserver reproducibility (30). Furthermore, cytokeratin IHC has been recommended by the International Tumor Bud Consensus Conference to help TB quantification in challenging cases such as tumor regions with strong inflammation or ruptured glands (2). In the current study, tumor budding was identified by the use of cytokeratin immunofluorescence and quantified by automated and standardized image analysis across the entire TBROI. Previous studies have highlighted the difficulty in the qualitative assessment of TBs in circumstances of gland fragmentation, stromal cells mimicking TBs, inflammation (45), microvesicles, and membrane fragments that might be cytokeratin-positive (30). In order to address these concerns, we ensured the presence of whole nuclei (n = 1–4) for the TB identification and only counted the isolated clusters of cells as buds. Quantifying TBs at the whole TBROI has the advantage of being representative of the whole IM rather than a manually selected region of interest. However, our methodology has the potential to “dilute” the significance of TBs represented by hotspots due to the tumor's innate heterogeneity. Like previous semiquantitative studies (6, 41), our quantification methodology reported TB number and density at the TBROI to be significantly correlated to poor patient survival.

The IMMUNOSCORE®*, proposed by Galon and colleagues (9, 46, 47), has been validated in a large international colon cancer study (14). It is not only a promising prognostic index (9, 11), but it has been shown to be superior to TNM staging when assessing colorectal cancer prognosis (46, 48). This score allows for the quantification of the in situ lymphocytic infiltration by the chromogenic labeling of two sequential tissue sections: one with CD3⁺ T cells and one with CD8⁺ T cells and quantified in two regions (IM and CT). The LIS created in this study, quantified lymphocytic infiltration from multiplexed immunofluorescence across a single tissue section. Our LIS results compared favorably with those of Galon and colleagues (9), showing that a high density of TILs was correlated with better survival outcomes in stage II colorectal cancer.

In order to reflect the dynamic processes present at the tumor front, T lymphocytes in close (50 and 100 μm) proximity of TBs in a two-dimensional plot were quantified. Our results revealed that a high number of T lymphocytes near TBs were associated with better survival. This supports the theory that TBs represent a heterogeneous subpopulation of cells and, although hypothesized to be more invasive than well-differentiated tumor glands, retain the potential to be recognized by immune-surveillance mechanisms. This process has previously been described as “nipping in the bud” (42). This theory was also supported by a study on the microenvironment of TBs by Lugli and colleagues (49), reporting high lymphocyte to TB ratio as a good prognostic factor in patients with colorectal cancer, and by Lang-Schwarz and colleagues (50), reporting the combination of both TILs and tumor budding predicting survival in colorectal cancer. Due to the high complexity and heterogeneity within the whole tissue slide, spatial mapping of specific features by eye is difficult, subjective, and time consuming. Selected regions of interest to study this interplay between different features and cells might not represent the spatial relationships present in the whole tissue. Here, we propose a rapid and completely unbiased method capturing every single event at the single-cell resolution and therefore better reporting on the interactions among TBs and TILs seen at the whole IM. Furthermore, such methodology can reveal specific patterns that may advance our understanding of the interactions occurring between cancer cells and their microenvironment.

Multiple features reported from the clinicopathologic report and by automated image analysis demonstrated prognostic significance. In order to reduce the dimensionality of the data and identify the independently significant features, a Cox regression model with LASSO regularization was applied to all 32 reported features in this study. This reported a reduced and independently significant set of eight features: four from the clinicopathologic report and four from image analysis. This reduced feature set was input for a RF machine learning model, which consisted of 500 trees. The advantage of feature selection through RF is the thorough testing of all significant features across multiple models while building in validation. This machine learning application reported only four features with a high mean decrease Gini coefficient, all of which were exported from image analysis. In order to avoid overfitting of our model on the training cohort, we included multiple validation steps within the workflow. Through this, we demonstrate that a model in which patients with a combined low CD3⁺ T-cell density in WTS, low mean number of TILs within 50 μm proximity to TBs, and high TB number conferred a survival disadvantage compared with any other combination. When a multivariate forward stepwise Cox regression including pT stage, TB number, LIS, and TBISI was applied, the only feature reported to significantly add value to the model was the TBISI, whereas pT stage, TB number, and LIS were not included in the equation. The new prognostic model of TBISI was then successfully validated across two independent cohorts, the second of which included patients with different genetic background, tumor genetic features, and gut microbiota, compared with the first cohort. This model and its cutoff points did not need to be altered for the Japanese validation cohort despite the fact that their tissue was sectioned longitudinally compared with the other cohorts that were cut at a cross-section. pT stage in the Japanese cohort, unlike both Edinburgh ones, was not significantly correlated with DSS in univariate analysis. This reflects the heterogeneous nature of stage II colorectal cancer disease and further highlights the need for additional methodologies, such as the TBISI reported here, to more accurately stratify these patients for both prognostic reasons and potentially identifying patients who may benefit from adjuvant therapy.

In conclusion, this study reports an automated image analysis–based workflow with the ability to quantify the tumor-infiltrating immune cells; total T cells (CD3⁺) and cytotoxic T cells (CD8⁺), and TBs at the invasive front, across a single tissue section. Additionally, through a machine learning approach, our results show evidence that the spatial association of lymphocytes and TBs has a high prognostic significance in stage II colorectal cancer patients. This combinational prognostic model demonstrated augmented accuracy and precision over that gained by either characterizing the lymphocytic infiltration or tumor budding in isolation. Furthermore, the final model contained no clinical parameters, thus demonstrating the benefit of automated profiling of the TME to provide a more precise and accurate prognosis. The methodology completely removes human subjectivity and negates the need for experienced staff to spend valuable time quantifying complex features across large areas of a whole slide section. This research can serve as a base for future studies on the prognostic significance of the interplay between different cell types within the patient's heterogeneous and heterotypic TME.

No potential conflicts of interest were disclosed.

Conception and design: I.P. Nearchou, P.D. Caie

Development of methodology: I.P. Nearchou, P.D. Caie

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): I.P. Nearchou, C.G. Gavriel, H. Ueno, D.J. Harrison, P.D. Caie

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): I.P. Nearchou

Writing, review, and/or revision of the manuscript: I.P. Nearchou, K. Lillard, H. Ueno, D.J. Harrison, P.D. Caie

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): I.P. Nearchou, K. Lillard

Study supervision: D.J. Harrison, P.D. Caie

The authors would like to acknowledge Frances Rae from Tissue Governance NHS Lothian and Dr. Yoshiki Kajiwara from the National Defense Medical College for securing ethics and patient follow-up data as well as In Hwa Um and Mustafa Elshani from the School of Medicine, University of St. Andrews, for cutting the tissue sections. Funding for the study was provided by Medical Research Scotland, and Indica Labs, Inc., provided in-kind resource.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.