Prognostic significance of myeloid immune cells and their spatial distribution in the colorectal cancer microenvironment.

BACKGROUND: Myeloid cells represent an abundant yet heterogeneous cell population in the colorectal cancer microenvironment, and their roles remain poorly understood.
METHODS: We used multiplexed immunofluorescence combined with digital image analysis to identify CD14+ monocytic and CD15+ granulocytic cells and to evaluate their maturity (HLA-DR and CD33), immunosuppressive potential (ARG1) and proximity to cytokeratin (KRT)-positive tumor cells in 913 colorectal carcinomas. Using covariate data of 4465 incident colorectal cancers in two prospective cohort studies, the inverse probability weighting method was used with multivariable-adjusted Cox proportional hazards models to assess cancer-specific mortality according to ordinal quartiles (Q1-Q4) of myeloid cell densities. Immune cell-tumor cell proximity was measured with the nearest neighbor method and the G-cross function, which determines the likelihood of any tumor cell having at least one immune cell of the specified type within a certain radius.
RESULTS: Higher intraepithelial (P trend=0.0002; HR for Q4 (vs Q1), 0.48, 95% CI 0.31 to 0.76) and stromal (P trend <0.0001; HR for Q4 (vs Q1), 0.42, 95% CI 0.29 to 0.63) densities of CD14+HLA-DR+ cells were associated with lower colorectal cancer-specific mortality while, conversely, higher intraepithelial densities of CD14+HLA-DR- cells were associated with higher colorectal cancer-specific mortality (P trend=0.0003; HR for Q4 (vs Q1), 1.78, 95% CI 1.25 to 2.55). Spatial analyses indicated that CD15+ cells were located closer to tumor cells than CD14+ cells, and CD14+HLA-DR+ cells were closer to tumor than CD14+HLA-DR- cells (p<0.0001). The G-cross proximity measurement, evaluating the difference in the likelihood of any tumor cell being colocated with at least one CD14+HLA-DR+ cell versus CD14+HLA-DR- cell within a 20 µm radius, was associated with lower colorectal cancer-specific mortality (P trend <0.0001; HR for Q4 (vs Q1), 0.37, 95% CI 0.24 to 0.57).
CONCLUSIONS: Myeloid cell populations occur in spatially distinct distributions and exhibit divergent, subset-specific prognostic significance in colorectal cancer, with mature CD14+HLA-DR+ and immature CD14+HLA-DR- monocytic phenotypes most notably showing opposite associations. These results highlight the prognostic utility of multimarker evaluation of myeloid cell infiltrates and reveal a previously unrecognized degree of spatial organization for myeloid cells in the immune microenvironment. © Author(s) (or their employer(s)) 2021. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

Background

Colorectal cancer is the second most common cause of cancer deaths in the world.1 A substantial fraction of cells present in the colorectal cancer microenvironment are derived from myeloid progenitor cells in the bone marrow.2 These cells are heterogenous, can exhibit functional and phenotypic plasticity and may exert protumorigenic and antitumorigenic effects.2 3 Analyses incorporating this diversity may advance the understanding of colorectal cancer biology and facilitate the development of improved tumor-immune biomarkers and cancer therapies.

In the late 1990s and early 2000s, it was established that cancer might cause aberrant, sustained myelopoiesis resulting in the accumulation of immature myeloid cells in the circulation and peripheral tissues4–6; such cells were subsequently termed myeloid-derived suppressor cells (MDSCs) in recognition of their high potential to suppress immune responses.7 Myeloid cells comprise CD14+ monocytic cells (such as monocytes and macrophages) and CD15+ granulocytic cells (such as neutrophils and eosinophils), and accordingly, the major human MDSC subgroups include CD14+ monocytic MDSCs and CD15+ polymorphonuclear MDSCs.6

Accurate identification of specific myeloid cell subsets requires a combination of multiple markers, which has been possible with multiparameter flow cytometry.6 Recent developments in multiplex immunofluorescence have enabled simultaneous detection of multiple antigens in situ with single cell resolution.8 However, these methods have not yet been used to analyze immune cells in the colorectal cancer microenvironment, where most such analyses have been conducted using traditional immunohistochemistry.2 9

In this study, we used multiplex immunofluorescence combined with digital image analysis and machine learning to identify CD14+ monocytic and CD15+ granulocytic myeloid cells and to evaluate their maturity (HLA-DR and CD33), suppressive potential (arginase 1 (ARG1)),10 and proximity to cytokeratin (KRT)-positive tumor cells in 913 colorectal cancers in two US nationwide prospective cohort studies. As our primary aim, using 4465 incident colorectal cancer cases and the inverse probability weighting (IPW) method to adjust for selection bias, we analyzed the prognostic significance of different myeloid cell populations within tumor tissue. As our secondary aim, we characterized the spatial patterns of different myeloid infiltrates and their prognostic significance, as well as the relationships of myeloid infiltrates with tumor molecular features.

Methods

Study population

The study was based on two prospective cohort studies, the Nurses’ Health Study (NHS; 121 701 women followed since 1976) and the Health Professionals Follow-up Study (HPFS; 51 529 men followed since 1986), where 4465 incident colorectal cancer cases had occurred and were documented until 2012. We included both colon and rectal carcinoma cases, based on the colorectal continuum model,11 and analyzed myeloid cell densities in 913 cases, based on the availability of adequate tissue specimens (table 1, online supplemental figure S1). All participants gave informed consent for the study.

Table 1

Clinical, pathological and molecular characteristics of colorectal cancer cases according to CD14+ cell and CD15+ cell densities

Characteristic*	All cases(N=913)	Overall CD14+ cell density				Overall CD15+ cell density
		Q1 (lowest)(N=228)	Q2(N=228)	Q3(N=229)	Q4 (highest)(N=228)	Q1 (lowest)(N=228)	Q2(N=228)	Q3(N=229)	Q4 (highest)(N=228)
Sex, n (%)
Female (NHS)	505 (55)	128 (56)	131 (57)	121 (53)	125 (55)	122 (54)	119 (52)	135 (59)	129 (57)
Male (HPFS)	408 (45)	100 (44)	97 (43)	108 (47)	103 (45)	106 (46)	109 (48)	94 (41)	99 (43)
Mean age±SD (years)	68.8±8.9	68.2±8.8	69.0±8.7	68.6±9.1	69.7±8.9	67.7±9.3	69.2±8.7	68.9±8.9	69.6±8.7
Year of diagnosis, n (%)
1995 or before	301 (33)	85 (37)	72 (32)	76 (33)	68 (30)	80 (35)	74 (32)	85 (37)	62 (27)
1996–2000	299 (33)	72 (32)	76 (33)	76 (33)	75 (33)	75 (33)	79 (35)	72 (31)	73 (32)
2001–2008	313 (34)	71 (31)	80 (35)	77 (34)	85 (37)	73 (32)	75 (33)	72 (31)	93 (41)
Family history of colorectal cancer in first-degree relative(s), n (%)
Absent	720 (79)	181 (80)	183 (81)	185 (81)	171 (76)	174 (77)	184 (81)	178 (78)	184 (81)
Present	189 (21)	46 (20)	44 (19)	44 (19)	55 (24)	51 (23)	43 (19)	51 (22)	44 (19)
Tumor location, n (%)
Cecum	159 (17)	33 (15)	31 (14)	46 (20)	49 (22)	38 (17)	32 (14)	37 (16)	52 (23)
Ascending to transverse colon	298 (33)	78 (34)	75 (33)	71 (31)	74 (33)	64 (28)	79 (35)	72 (32)	83 (36)
Descending to sigmoid colon	272 (30)	71 (31)	69 (30)	69 (30)	63 (28)	82 (36)	64 (28)	75 (33)	51 (22)
Rectum	178 (20)	45 (20)	53 (23)	41 (18)	41 (18)	42 (19)	52 (23)	44 (19)	42 (18)
AJCC disease stage, n (%)
I	193 (23)	44 (20)	47 (22)	52 (24)	50 (24)	40 (19)	45 (21)	55 (26)	53 (24)
II	283 (33)	67 (31)	65 (31)	69 (32)	82 (39)	66 (31)	71 (34)	73 (35)	73 (34)
III	245 (29)	66 (30)	64 (30)	59 (28)	56 (27)	64 (30)	62 (30)	54 (26)	65 (30)
IV	130 (15)	41 (19)	36 (17)	33 (15)	20 (9.6)	43 (20)	32 (15)	29 (14)	26 (12)
Tumor differentiation, n (%)
Well to moderate	824 (90)	212 (93)	215 (94)	204 (89)	193 (85)	209 (92)	212 (93)	211 (92)	192 (84)
Poor	88 (9.7)	15 (6.6)	13 (5.7)	25 (11)	35 (15)	18 (7.9)	16 (7.0)	18 (7.9)	36 (16)
MSI status, n (%)
Non-MSI-high	735 (83)	196 (87)	190 (87)	180 (81)	169 (75)	209 (94)	194 (86)	179 (82)	153 (69)
MSI-high	154 (17)	29 (13)	28 (13)	42 (19)	55 (25)	14 (6.3)	31 (14)	39 (18)	70 (31)
CIMP status, n (%)
Low/negative	686 (81)	187 (87)	188 (87)	160 (77)	151 (73)	188 (90)	181 (85)	170 (81)	147 (70)
High	159 (19)	29 (13)	27 (13)	48 (23)	55 (27)	22 (10)	32 (15)	41 (19)	64 (30)
Mean LINE-1 methylation level±SD (%)	62.5±9.5	61.5±9.6	61.6±9.7	62.5±9.5	64.5±9.1	61.1±9.8	61.5±9.5	62.8±9.6	64.8±8.9
KRAS mutation, n (%)
Wild-type	530 (60)	121 (54)	134 (61)	127 (57)	148 (67)	129 (58)	128 (57)	135 (62)	138 (63)
Mutant	356 (40)	102 (46)	87 (39)	94 (43)	73 (33)	93 (42)	98 (43)	84 (38)	81 (37)
BRAF mutation, n (%)
Wild-type	757 (85)	192 (86)	196 (87)	181 (81)	188 (85)	194 (86)	198 (89)	192 (86)	173 (78)
Mutant	136 (15)	32 (14)	28 (13)	43 (19)	33 (15)	31 (14)	24 (11)	31 (14)	50 (22)
PIK3CA mutation, n (%)
Wild-type	698 (84)	179 (84)	173 (85)	182 (86)	164 (80)	187 (89)	171 (81)	175 (86)	165 (79)
Mutant	136 (16)	35 (16)	30 (15)	29 (14)	42 (20)	24 (11)	41 (19)	28 (14)	43 (21)
Neoantigen load, n (%)
Q1 (lowest)	103 (25)	30 (29)	22 (21)	25 (24)	26 (25)	28 (27)	25 (24)	29 (28)	21 (20)
Q2	103 (25)	29 (28)	33 (32)	23 (22)	18 (17)	30 (29)	25 (24)	24 (23)	24 (23)
Q3	103 (25)	23 (22)	22 (21)	32 (31)	26 (25)	27 (26)	27 (26)	22 (21)	27 (26)
Q4 (highest)	103 (25)	21 (20)	26 (25)	23 (22)	33 (32)	18 (17)	26 (25)	28 (27)	31 (30)

*Percentage indicates the proportion of patients with a specific clinical, pathologic or molecular characteristic among all patients or in strata of myeloid cell density.

AJCC, American Joint Committee on Cancer; CIMP, CpG island methylator phenotype; HPFS, Health Professionals Follow-up Study; LINE-1, long-interspersed nucleotide element-1; MSI, microsatellite instability; NHS, Nurses’ Health Study.;

Blinded to exposure data, study physicians reviewed the medical records related to colorectal cancer and collected clinical information such as the American Joint Committee on Cancer tumor, node, metastases stage and tumor location. Deaths and unreported lethal colorectal cancer cases were identified using the National Death Index, and cause of death was also recorded. Survival time was defined as the period from the date of colorectal cancer diagnosis to death or the end of follow-up (January 1, 2016 for HPFS; June 1, 2016 for NHS). For analyses of colorectal cancer specific mortality, deaths as a result of other causes were censored.

Formalin-fixed paraffin-embedded (FFPE) tissue blocks were collected, and DNA was extracted using the QIAmp DNA Mini Kit (Qiagen, Valencia, California, USA).12 Microsatellite instability (MSI) status was analyzed with 10 microsatellite markers (BAT25, BAT26, BAT40, D2S123, D5S346, D17S250, D18S55, D18S56, D18S67, and D18S487), as previously described.13 Methylation analyses for eight promoters specific for the CpG island methylator phenotype (CIMP; CACNA1G, CDKN2A, CRABP1, IGF2, MLH1, NEUROG1, RUNX3, and SOCS1) were conducted using bisulfite-treated DNA and real-time PCR,13 and long interspersed nucleotide element-1 (LINE-1) methylation level was assessed with pyrosequencing.14 Pyrosequencing was performed to evaluate BRAF (codon 600), KRAS (codons 12, 13, 61 and 146) and PIK3CA (exons 9 and 20) mutation status.15 Neoantigen load of 412 tumors was estimated using a neoantigen prediction pipeline for somatic mutations based on whole-exome sequencing and identifying peptides binding to human leukocyte antigen (MHC) molecules with high affinity.16 One pathologist (SO) reviewed H&E-stained tissue sections to confirm the diagnoses and evaluate tumor differentiation (well/moderate vs poor) and lymphocytic reaction patterns (tumor-infiltrating lymphocytes; intratumoral periglandular reaction; peritumoral reaction; and Crohn’s-like lymphoid reaction).17 Tissue microarrays (TMAs) containing 1–4 (median 2) 0.6 mm diameter cores from representative areas of each tumor were constructed.18 Densities of granulocytic cells (neutrophils and eosinophils) from 900 tumors were evaluated from H&E-stained TMA sections using a computer-assisted method, as previously described.19

Multiplex immunofluorescence

We built a multiplex immunofluorescence assay, incorporating tyramide signal amplification for high antigen sensitivity,20 to evaluate myeloid cells in the colorectal cancer microenvironment. First, we screened potential markers and antibody clones based on published literature.6 21–27 Candidate antibody testing was then performed on colorectal cancer tissue and lymphoid controls via chromogenic immunohistochemistry, and antibody clones with expression patterns corresponding to their biologically expected distribution (at both the tissue level and subcellular level) were selected for further testing via immunofluorescence. Antibody clones with the highest signal-to-noise ratio were then compiled into a multiplex immunofluorescence panel containing a monocytic lineage marker (CD14), a granulocytic lineage marker (CD15), a macrophage maturity marker (HLA-DR), a myeloid cell immaturity marker (CD33), a suppressive marker (ARG1), and a tumor epithelial cell marker (KRT, keratins), following standardized protein nomenclature recommended by a panel of experts28 (online supplemental figure S2). Antibody and fluorophore concentrations and antibody–fluorophore pairing were then iteratively optimized, and similar staining patterns between chromogenic immunohistochemistry and multiplex immunofluorescence were confirmed (online supplemental figure S3). The final immunofluorescence protocol was automated with a Leica Bond RX Research Stainer (Leica Biosystems, Buffalo, Illinois, USA) (online supplemental figure S2).

Stained slides were imaged using the Vectra 3.0 multispectral imaging system (Akoya Biosciences, Hopkinton, Massachusetts, USA) equipped with a 20× objective.29 30 TMA core images were reviewed, and images containing more than 10% tumor across their total area were processed via pathologist-supervised machine learning algorithms using the inForm software package (V.2.4.1; Akoya Biosciences) to classify tissue compartments into three categories (tumor epithelium, stroma and other) and to detect and classify cells into four main categories (CD14+, CD15+, tumor cell and other) (figure 1A–D).29 31

Figure 1

Quantitative, multiplexed analysis of myeloid immune cell populations and their spatial localization. (A) Multiplex immunofluorescence image. Scale bar is 100 µm. The images were analyzed using supervised machine learning algorithms including tissue category classification (B) cell segmentation (C), and cell type classification (D). (E) Distances from immune cells to the closest tumor cell (NNDMyeloid cell:Tumor) were calculated. (F) NNDMyeloid cell:Tumor according to the main myeloid cell categories. P values are for the Wilcoxon rank sum test. (G) Marker signal intensities according to NNDMyeloid cell:Tumor. ****P<0.0001, ***P<0.001, NS: P>0.005.

Single cell-level data were exported from inForm and further processed with the R statistical programming language (V.3.5.2; R Foundation for Statistical Computing, Vienna, Austria). Marker intensities were scaled by calculating Z-scores within each TMA, in order to correct for potential staining variability, and myeloid cell subpopulations were defined based on scaled mean intensities of HLA-DR (membrane), CD33 (membrane), and ARG1 (cytoplasm) expression in specific cellular compartments of CD14+ monocytic lineage cells and CD15+ granulocytic lineage cells. We conducted spatial point pattern analysis using the spatstat R package.32 For each myeloid cell, we computed the distance to the nearest tumor cells (NNDMyeloid cell:Tumor) (figure 1E). For visualization, fluorophore intensities of individual cells across all cores (scaled within the cell population-of-interest) as a function of NNDMyeloid cell:Tumor were plotted with the ggplot2 package using generalized additive model smoothing (formula=y ~ s(x)). We used the G-cross function (GTumor:Myeloid cell(r)) to estimate the probability at least one immune cell being located within a radius r of any tumor cell. For each sample, we calculated the empirical GTumor:Myeloid cell(r) distribution functions based on observed nearest neighbor distances, applying edge correction via the Kaplan-Meier method. For survival analysis, we evaluated function values at a radius of 20 µm (GTumor:Myeloid cell(20 µm)). This radius was preselected prior to analyses in order to identify immune cell populations likely capable of effective, direct, cell-to-cell interaction with tumor cells, consistent with prior reports in multiple tumor types, including colorectal cancer.19 33

Statistical analysis

Statistical analyses were performed using SAS software (V.9.4, SAS Institute, Cary, North Carolina, USA). Our primary hypothesis tested the relationship between each myeloid cell density variable and colorectal cancer-specific survival using multivariable-adjusted Cox proportional hazards regression analysis. We used the stringent two-sided α level of 0.005 to judge statistical significance.34 All other analyses were secondary, and we interpreted the results cautiously.

We evaluated relationships between myeloid cell densities and clinicopathologic features using the χ2 test and Spearman’s rank correlation (rho) test, as appropriate. We used the Kaplan-Meier method to estimate survival functions and the log-rank test to compare survival distributions. To adjust for selection bias due to tissue availability, the IPW method35–37 was used with multivariable-adjusted Cox proportional hazards regression models that assessed the HRs for colorectal cancer-specific mortality according to ordinal myeloid cell density categories (Q1–Q4). Further details of the statistical methods used in the survival analyses are shown in online supplemental table S1.

Results

Characterization of myeloid infiltrates in colorectal cancer

We evaluated myeloid immune infiltrates in 913 colorectal cancer cases (table 1). Immunofluorescence images were processed with pathologist-supervised machine learning algorithms (figure 1A–D), generating data for 2 541 019 tumor cells, 6 63 854 CD14+ monocytic cells, and 131 470 CD15+ granulocytic cells. Fluorophore signal intensity distributions were comparable across the TMAs (online supplemental figure S4). Core-to-core correlation for the densities of both overall CD14+ cells (online supplemental figure S5) and CD15+ cells (online supplemental figure S6) were reasonably high (Spearman rho=0.50–0.65). CD15+ cell density was moderately correlated with granulocytic cell density as evaluated by digital image analysis of H&E-stained sections (Spearman rho=0.44).

Overall, CD14+ monocytic cells (median density 994 cells/mm2, IQR 609–1580 cells/mm2) were more common than CD15+ granulocytic cells (median density 97 cells/mm2, IQR 28–256 cells/mm2). Densities of both CD14+ cells and CD15+ cells were higher in tumor stroma than within tumor intraepithelial regions (median stromal:intraepithelial density ratios 12.5 and 5.0, respectively). The majority (75%) of CD14+ monocytic cells expressed the antigen presentation molecule HLA-DR, indicating a mature macrophage phenotype (online supplemental figure S7), although expression intensity was variable. In contrast, expression of myeloid progenitor marker CD33 was relatively rare (<10%) in both CD14+ and CD15+ cell lineages (online supplemental figure S7). Expression of immunosuppressive enzyme ARG1 was common in CD15+ cells (>90%) yet infrequent in CD14+ cells (<10%) (online supplemental figure S7). Intraepithelial and stromal densities of CD14+ cells (Spearman rho=0.55) (or those of CD15+ cells (Spearman rho=0.81)) were positively correlated with each other (online supplemental figure S8), while densities of CD14+ cells and CD15+ cells showed lower correlation (intraepithelial: Spearman rho=0.33; stromal: Spearman rho=0.22) (online supplemental figure S8).

We calculated the distance from each myeloid cell to the closest tumor cell (figure 1E). We found that CD15+ cells were located, on average, closer to tumor cells than CD14+ cells (figure 1F–G). Within the CD14+ monocytic population, mature HLA-DR+ cells were located closer to tumor cells, while immature CD14+CD33+ cells were less abundant near tumor cells (figure 1F–G). Within the CD15+ granulocytic population, ARG1 expression was higher in cells located close to tumor cells (figure 1F–G). Additional example images depicting these myeloid cell spatial distributions are shown in online supplemental figure S9.

Survival analyses

As our primary hypothesis, we evaluated the prognostic significance of myeloid immune cell densities using a cohort of 4465 incident colorectal carcinoma cases that included the 913 cases with available myeloid cell data and the IPW method to adjust for selection bias due to tissue availability. During the median follow-up time of 16.3 years (IQR 12.8–20.1 years) for censored cases, there were 613 all-cause deaths, including 283 colorectal cancer-specific deaths.

We found that higher stromal densities of CD14+ cells (P trend <0.0001; HR for Q4 (vs Q1) of 0.47, 95% CI 0.31 to 0.69) and CD15+ cells (P trend=0.0036; HR for Q4 (vs Q1) of 0.61, 95% CI 0.42 to 0.89) were associated with longer cancer-specific survival independent of potential confounding factors such as disease stage and MSI status (table 2, figure 2 and online supplemental figure S10). In contrast, intraepithelial CD14+ or CD15+ cell densities were not statistically significant (P trend >0.01 in multivariable models at the α level of 0.005), highlighting the value of studying immune cell densities in defined tissue regions.

Table 2

Densities of CD14+ cells and CD15+ cells in tumor intraepithelial, stromal, and overall regions and patient survival with inverse probability weighting (IPW)

	No. of cases	Colorectal cancer-specific survival			Overall survival
		No. of events	UnivariableHR (95% CI)*	MultivariableHR (95% CI)*†	No. of events	UnivariableHR (95% CI)*	MultivariableHR (95% CI)*†
Tumor intraepithelial region
CD14+ cell density
Q1	228	88	1 (referent)	1 (referent)	158	1 (referent)	1 (referent)
Q2	228	64	0.67 (0.47 to 0.95)	0.79 (0.56 to 1.11)	153	0.91 (0.68 to 1.20)	1.01 (0.76 to 1.33)
Q3	229	67	0.75 (0.53 to 1.05)	0.74 (0.51 to 1.06)	151	0.81 (0.60 to 1.09)	0.79 (0.58 to 1.08)
Q4	228	64	0.71 (0.50 to 1.01)	0.82 (0.55 to 1.22)	151	0.89 (0.66 to 1.18)	0.91 (0.64 to 1.28)
P trend‡			0.086	0.22		0.30	0.29
CD15+ cell density
Q1	228	86	1 (referent)	1 (referent)	157	1 (referent)	1 (referent)
Q2	228	67	0.76 (0.54 to 1.08)	0.97 (0.69 to 1.35)	156	0.97 (0.73 to 1.30)	1.09 (0.82 to 1.45)
Q3	229	77	0.82 (0.58 to 1.15)	0.94 (0.66 to 1.35)	162	0.99 (0.74 to 1.33)	1.04 (0.76 to 1.42)
Q4	228	53	0.56 (0.39 to 0.82)	0.58 (0.39 to 0.86)	138	0.77 (0.57 to 1.05)	0.76 (0.55 to 1.06)
P trend‡			0.0061	0.011		0.14	0.12
Tumor stromal region
CD14+ cell density
Q1	228	93	1 (referent)	1 (referent)	169	1 (referent)	1 (referent)
Q2	228	86	0.90 (0.66 to 1.24)	1.07 (0.79 to 1.46)	158	0.96 (0.73 to 1.27)	1.09 (0.83 to 1.44)
Q3	229	59	0.61 (0.42 to 0.86)	0.68 (0.47 to 0.97)	145	0.73 (0.55 to 0.98)	0.77 (0.57 to 1.04)
Q4	228	45	0.38 (0.26 to 0.56)	0.47 (0.31 to 0.69)	141	0.55 (0.41 to 0.73)	0.53 (0.39 to 0.73)
P trend‡			<0.0001	<0.0001		<0.0001	<0.0001
CD15+ cell density
Q1	228	95	1 (referent)	1 (referent)	160	1 (referent)	1 (referent)
Q2	228	71	0.67 (0.48 to 0.93)	0.72 (0.53 to 1.00)	163	0.79 (0.60 to 1.05)	0.77 (0.58 to 1.02)
Q3	229	67	0.50 (0.35 to 0.70)	0.60 (0.42 to 0.86)	153	0.65 (0.48 to 0.87)	0.67 (0.49 to 0.91)
Q4	228	50	0.51 (0.35 to 0.75)	0.61 (0.42 to 0.89)	137	0.67 (0.50 to 0.90)	0.71 (0.52 to 0.97)
P trend‡			0.0001	0.0036		0.0039	0.016
Overall region
CD14+ cell density
Q1	228	83	1 (referent)	1 (referent)	161	1 (referent)	1 (referent)
Q2	228	73	0.83 (0.59 to 1.18)	0.87 (0.62 to 1.24)	155	0.96 (0.72 to 1.28)	0.99 (0.74 to 1.32)
Q3	229	74	0.88 (0.63 to 1.24)	0.94 (0.66 to 1.33)	156	0.99 (0.75 to 1.33)	0.95 (0.70 to 1.29)
Q4	228	53	0.57 (0.39 to 0.83)	0.68 (0.46 to 1.01)	141	0.74 (0.55 to 1.00)	0.70 (0.51 to 0.96)
P trend‡			0.0074	0.097		0.075	0.034
CD15+ cell density
Q1	228	88	1 (referent)	1 (referent)	157	1 (referent)	1 (referent)
Q2	228	73	0.78 (0.55 to 1.09)	0.86 (0.61 to 1.20)	162	0.97 (0.73 to 1.29)	0.90 (0.68 to 1.19)
Q3	229	67	0.59 (0.41 to 0.84)	0.66 (0.46 to 0.95)	158	0.72 (0.53 to 0.97)	0.70 (0.51 to 0.95)
Q4	228	55	0.63 (0.44 to 0.91)	0.80 (0.54 to 1.17)	136	0.80 (0.59 to 1.08)	0.84 (0.61 to 1.15)
P trend‡			0.0048	0.092		0.043	0.099

*IPW was applied to reduce bias due to the availability of tumor tissue after cancer diagnosis (see ‘Statistical analysis’ subsection for details).

†The multivariable Cox regression model initially included sex, age, year of diagnosis, family history of colorectal cancer, tumor location, tumor differentiation, disease stage, microsatellite instability, CpG island methylator phenotype, KRAS, BRAF, and PIK3CA mutations, and long-interspersed nucleotide element-1 methylation level. A backward elimination with a threshold P of 0.05 was used to select variables for the final models.

‡The P trend value was calculated across the ordinal quartiles of each myeloid cell density within tumor epithelial, stromal, and overall regions in the IPW-adjusted Cox regression model.

Figure 2

Inverse probability weighting-adjusted Kaplan-Meier survival curves for colorectal cancer-specific survival according to ordinal quartile categories (Q1–Q4) of intraepithelial (A–C) and stromal (D–F) densities of CD14+ cells classified by HLA-DR expression.

Having identified an association between CD14+ monocytic cell density and longer survival, we sought to further classify these cells according to maturity using HLA-DR expression. This revealed a significant association between prognosis and degree of maturity (table 3, figure 2). Higher densities of both intraepithelial (P trend=0.0002; HR for Q4 (vs Q1) of 0.48, 95% CI 0.31 to 0.76) and stromal (P trend <0.0001; HR for Q4 (vs Q1) of 0.42, 95% CI 0.29 to 0.63) CD14+HLA-DR+ mature monocytic cells were associated with lower cancer-specific mortality while, conversely, higher intraepithelial densities of CD14+HLA-DR− immature monocytic cells were associated with higher cancer-specific mortality (P trend=0.0003; HR for Q4 (vs Q1) of 1.78, 95% CI 1.25 to 2.55). Given that MSI status strongly influences the antitumor immune response,38 we further examined whether the prognostic association of CD14+HLA-DR+ or CD14+HLA-DR− myeloid cell densities might differ by MSI status and found that the difference was not statistically significant at an α level of 0.005 (P interaction >0.03) (online supplemental table S3).

Table 3

Densities of CD14+ cells classified according to HLA-DR expression in tumor intraepithelial and stromal regions and patient survival with inverse probability weighting (IPW)

	No. of cases	Colorectal cancer-specific survival			Overall survival
		No. of events	UnivariableHR (95% CI)*	MultivariableHR (95% CI)*†	No. of events	UnivariableHR (95% CI)*	MultivariableHR (95% CI)*†
Tumor intraepithelial region
CD14+HLA-DR+ cell density
Q1	228	94	1 (referent)	1 (referent)	158	1 (referent)	1 (referent)
Q2	228	76	0.88 (0.63 to 1.21)	0.90 (0.66 to 1.23)	157	1.08 (0.82 to 1.43)	1.06 (0.81 to 1.40)
Q3	229	65	0.66 (0.47 to 0.94)	0.64 (0.45 to 0.93)	161	0.82 (0.61 to 1.11)	0.76 (0.56 to 1.04)
Q4	228	48	0.44 (0.30 to 0.64)	0.48 (0.31 to 0.76)	137	0.64 (0.48 to 0.87)	0.61 (0.42 to 0.87)
P trend‡			<0.0001	0.0002		0.0010	0.0011
CD14+HLA-DR− cell density
Q1	228	59	1 (referent)	1 (referent)	144	1 (referent)	1 (referent)
Q2	228	56	1.04 (0.69 to 1.56)	1.00 (0.67 to 1.49)	150	1.11 (0.82 to 1.51)	1.07 (0.78 to 1.48)
Q3	229	68	1.39 (0.95 to 2.05)	1.42 (0.97 to 2.09)	148	1.30 (0.96 to 1.76)	1.29 (0.96 to 1.75)
Q4	228	100	2.17 (1.52 to 3.08)	1.78 (1.25 to 2.55)	171	1.84 (1.37 to 2.47)	1.59 (1.18 to 2.15)
P trend‡			<0.0001	0.0003		<0.0001	0.0010
Tumor stromal region
CD14+HLA-DR+ cell density
Q1	228	107	1 (referent)	1 (referent)	177	1 (referent)	1 (referent)
Q2	228	76	0.66 (0.48 to 0.91)	0.77 (0.56 to 1.06)	154	0.71 (0.53 to 0.93)	0.76 (0.57 to 1.00)
Q3	229	58	0.45 (0.32 to 0.65)	0.52 (0.36 to 0.75)	142	0.57 (0.43 to 0.77)	0.59 (0.44 to 0.80)
Q4	228	42	0.33 (0.23 to 0.49)	0.42 (0.29 to 0.63)	140	0.51 (0.38 to 0.67)	0.51 (0.38 to 0.69)
P trend‡			<0.0001	<0.0001		<0.0001	<0.0001
CD14+HLA-DR− cell density
Q1	228	69	1 (referent)	1 (referent)	152	1 (referent)	1 (referent)
Q2	228	56	0.87 (0.59 to 1.28)	1.01 (0.69 to 1.48)	144	1.01 (0.75 to 1.37)	1.07 (0.78 to 1.48)
Q3	229	82	1.34 (0.93 to 1.91)	1.33 (0.92 to 1.91)	161	1.35 (1.00 to 1.81)	1.35 (0.99 to 1.84)
Q4	228	76	1.10 (0.76 to 1.58)	1.24 (0.86 to 1.78)	156	1.11 (0.82 to 1.49)	1.18 (0.86 to 1.62)
P trend‡			0.24	0.12		0.24	0.15

*IPW was applied to reduce bias due to the availability of tumor tissue after cancer diagnosis (see ’Statistical analysis’ subsection for details).

†The multivariable Cox regression model initially included sex, age, year of diagnosis, family history of colorectal cancer, tumor location, tumor differentiation, disease stage, microsatellite instability, CpG island methylator phenotype, KRAS, BRAF, and PIK3CA mutations, and long-interspersed nucleotide element-1 methylation level. A backward elimination with a threshold P of 0.05 was used to select variables for the final models.

‡The P trend value was calculated across the ordinal quartiles of each myeloid cell density within tumor epithelial and stromal regions in the IPW-adjusted Cox regression model.

Densities of immature CD14+ cells or CD15+ cells expressing CD33 were not significantly associated with cancer-specific survival, irrespective of stromal or intraepithelial location. In contrast, the prognostic associations of CD14+CD33− cell or CD15+CD33− cell densities were comparable to those of overall CD14+ or CD15+ cell densities, respectively, as most CD14+ cells and CD15+ cells were CD33− (online supplemental table S4). CD14+ or CD15+ cells expressing immunosuppressive marker ARG1 were not significantly associated with mortality (online supplemental table S5).

Relationships with tumor and patient characteristics

In order to understand the factors that govern myeloid cell infiltration in CRC, we analyzed the relationships between myeloid cell densities and clinicopathologic features (online supplemental table S11). Higher densities of both intraepithelial CD14+ cells (including HLA-DR+ and HLA-DR− subpopulations) and CD15+ cells were associated with MSI high status and poor differentiation (p<0.0002). Higher intraepithelial CD14+ and CD14+HLA-DR+ cell densities were also associated with higher neoantigen load (p<0.0002). In stroma, higher densities of CD14+ and CD14+HLA-DR+ cells were associated with lower stage (p<0.0001). Additionally, higher stromal densities of CD14+ cells, CD14+HLA-DR+ cells and CD15+ cells (p<0.003) were associated with MSI-high status.

Considering the opposing survival associations for CD14+HLA-DR+ and CD14+HLA-DR− cells, we assessed whether these populations might also be differentially associated with lymphocytic reaction patterns (online supplemental figure S12). We found that higher intraepithelial and stromal densities of CD14+HLA-DR+ cells were associated with more intense lymphocytic reaction patterns (tumor-infiltrating lymphocytes, intratumoral periglandular reaction, peritumoral reaction, and Crohn’s-like lymphoid reaction) (all p<0.0001), while CD14+HLA-DR− cell density was not significantly associated with any lymphocytic reaction pattern (p>0.02 with α=0.005) (online supplemental figure S12).

Spatial analyses with G-cross function

Given our ability to precisely define the positions of individual tumor cells and myeloid cells, we sought to evaluate the prognostic significance of proximity between tumor cells and myeloid cells. We used the GTumor:Myeloid cell(r) function to evaluate the likelihood of any tumor cell in the sample having at least one myeloid cell (of the specified type) within a specified radius r. Higher function values at a given radius r indicate greater colocalization of the specified myeloid cell type with tumor cells (examples in figure 3A–B, online supplemental figure S13A–B).

Figure 3

Spatial analysis of CD14+HLA-DR+ and CD14+HLA-DR− myeloid immune infiltrates with the tumor:myeloid cell G-cross function (GTumor:Myeloid cell). (A and B) Example myeloid cell patterns and corresponding GTumor:Myeloid cell (r) plots, estimating the probability of any tumor cell having at least one neighboring myeloid cell of the specified type within an r µm radius. (C) Univariable (black) and multivariable (red) Cox proportional hazards regression models for cancer-specific survival according to GTumor:Myeloid cell (20 µm) ordinal quartile categories (C1–C4). The multivariable Cox regression models initially included sex, age, year of diagnosis, family history of colorectal cancer, tumor location, tumor differentiation, disease stage, microsatellite instability, CpG island methylator phenotype, KRAS, BRAF, and PIK3CA mutations, and long-interspersed nucleotide element-1 methylation level. A backward elimination with a threshold P of 0.05 was used to select variables for the final models.

We found that the values of functions evaluating overall CD14+ and CD15+ populations within a 20 µm radius of tumor cells (GTumor:CD14+(20 µm) and GTumor:CD15+(20 µm)) were not significantly associated with cancer-specific survival (P trend >0.1) (online supplemental figure S13). Conversely, the values of functions evaluating both mature CD14+HLA-DR+ (P trend=0.0004) and immature CD14+HLA-DR− subsets (P trend=0.0007) were significantly prognostic and showed opposite effects (figure 3). Accordingly, their difference (GTumor:CD14+HLADR+(20 µm)-GTumor:CD14+HLADR−(20 µm)), where high values represent a high likelihood of any tumor cell being colocated with at least one mature monocytic cell yet a low likelihood of colocation with an immature monocytic cell, exhibited an even stronger prognostic association (P trend <0.0001; HR for Q4 (vs Q1), 0.37, 95% CI 0.24 to 0.57). Given that this association could be confounded by cell density, we directly compared the prognostic power of mature and immature monocytic cell densities and the G-cross proximity measurement by including the proximity and density variables in one Cox regression model, with reciprocal adjustment (online supplemental table S6). This analysis indicated that both CD14+HLA-DR+ and CD14+HLA-DR− cell densities and the G-cross proximity measurement independently contributed to prognostic significance.

Discussion

We leveraged recent technical advances in multiplexed immunofluorescence8 to build a novel assay enabling in situ characterization of detailed myeloid cell phenotypes in colorectal cancer. In addition to monocytic (CD14) and granulocytic (CD15) lineage markers, our panel contained macrophage maturation marker HLA-DR (an MHC class II cell surface receptor that supports antigen presentation to T cells),39 myeloid marker CD33 (a transmembrane sialic acid binding receptor downregulated during myeloid cell maturation),40 immunosuppressive marker ARG1 (an enzyme catalyzing the hydrolysis of arginine to ornithine and urea),10 and pan-cytokeratin (KRT) as a tumor cell marker. This panel design enabled us to simultaneously examine myeloid cell lineage (monocytic vs granulocytic), maturity (HLA-DR and CD33), and suppressive potential (ARG1) in a manner not possible with single marker approaches. This combinatorial approach was required for accurate cell identification since most of the included markers are not specific for the myeloid cell populations under study. For example, HLA-DR can be expressed by B cells and dendritic cells in addition to monocytic lineage cells,39 41 while CD33 is expressed by mast cells.22

We applied our assay to 913 colorectal cancers in two US nationwide prospective cohort studies. We found that high CD14+ monocytic and CD15+ granulocytic cell densities in tumor stroma were associated with low cancer-specific mortality independent of potential confounding factors such as disease stage, MSI status, and BRAF, KRAS, and PIK3CA mutation status. Furthermore, we found that these populations were heterogenous, and finer-grained classification of CD14+ monocytic cells according to HLA-DR expression revealed opposing prognostic roles for mature HLA-DR+ and immature HLA-DR− subsets within intraepithelial regions, with high density of mature HLA-DR+ cells associating with lower cancer-specific mortality and high density of intraepithelial immature HLA-DR− cells associating with higher cancer-specific mortality.

Since the discovery of the expansion of immature myeloid cells in cancer,4 5 a multitude of studies have assessed the potential mechanisms involved in their tumor-promoting activity.6 42 Experimental evidence suggests that these cells can produce anti-inflammatory cytokines, including interleukin 1043 and transforming growth factor beta 1,44 thereby reducing T cell infiltration.45 Moreover, they may produce angiogenic mediators or growth factors promoting tumor cell proliferation.45 46 Conversely, the mechanism linking mature HLA-DR+ monocytic cells with better survival may be related to this cell type’s ability to efficiently present antigens to T lymphocytes,47 thereby activating antitumor immunity. Indeed, we found that high densities of CD14+HLA-DR+ cells, but not CD14+HLA-DR− cells, were associated with high lymphocytic reaction scores.

In addition to CD14+ monocytic cells, higher CD15+ granulocytic cell density in tumor stromal regions was associated with favorable prognosis. This is consistent with several previous studies that have evaluated the prognostic significance of CEACAM8+ (CD66b+)26 48 49 granulocytic cells in colorectal cancer and stands in contrast to other tumor types, such as hepatocellular carcinoma50 and renal cell carcinoma,51 where strong granulocytic infiltrates have been associated with worse prognosis. These tumor type-specific results may be attributable to the diversity and plasticity of granulocytes within the tumor microenvironment, as such cells can participate in both tumor-promoting functions (such as angiogenesis, extracellular matrix remodeling, and immunosuppression) and antitumor roles (such as direct killing of tumor cells).45 Our finding of the favorable prognostic value of CD15+ granulocytes suggests that the antitumorigenic granulocytic subpopulations may outnumber protumorigenic ones within the colorectal cancer microenvironment.

ARG1 is considered one of the most important immunoregulatory enzymes in the tumor microenvironment, contributing to suppression of the T cell response through arginine depletion.10 52 In murine tumor models, high ARG1 expression has commonly been observed in monocytic lineage cells whereas, in humans, strong expression is most frequently seen on polymorphonuclear cells instead of monocytic cells.10 53 Consistent with these prior human studies, we found that ARG1 positivity was common in CD15+ granulocytic cells but relatively rare in CD14+ monocytic cells in the colorectal cancer microenvironment in a large human population-based sample. Contrary to our exploratory hypothesis, higher densities of granulocytic or monocytic ARG1+ cells were not associated with adverse clinical outcome. This may be related to: (1) a lower functional relevance for ARG1 as compared with other immunosuppressive pathways in colorectal cancer (such as the IDO1 pathway and the PDCD1 (PD-1)-CD274 (PD-L1) pathway,52 54 which were not evaluated in this study), (2) potential production of proinflammatory mediators by ARG1+ cells, thereby activating the immune response, or (3) additional unidentified factors, all of which are relevant topics for future studies. Of note, murine and human immune systems may have significant differences55; many promising findings in mouse models have not translated into the treatment of human disease, supporting the clinical relevance of large-scale analyses of the tumor microenvironment using human specimens.

We discovered previously unappreciated differences in the infiltration patterns of various myeloid subsets. Both CD14+ monocytic cells and CD15+ granulocytic cells were preferably located in tumor stroma, although CD15+ cells were, on average, closer to tumor cells than CD14+ cells. Within the CD14+ population, mature CD14+HLA-DR+ cells (and CD14+CD33− cells) were located closer to tumor cells as compared with immature CD14+HLA-DR− cells (and CD14+CD33+ cells). This suggests that the maturity of CD14+ cells may influence their interaction with tumor cells, with the HLA-DR+ subset having closer contact with tumor cells as compared with the HLA-DR− subset. Alternatively, close contact with tumor cells may itself drive the maturation of CD14+ cells. We hypothesize that these mature monocytic cells may participate in tumor cell phagocytosis, considering their association with favorable prognosis. Notably, however, CD14+HLA-DR+ cells were associated with longer survival in both intraepithelial and stromal compartments, while the adverse prognostic effect of CD14+HLA-DR− cells was only seen in the intraepithelial compartment, suggesting that the significance of the immature HLA-DR− subset rather than the mature HLA-DR+ subset may be dependent on the close interaction with tumor cells. Modeling the proximity of tumor cells with CD14+HLA-DR+ and CD14+HLA-DR− cells with G-cross function, we found that the G-cross proximity measurement (evaluating the difference in the likelihood of any tumor cell in the sample having at least one CD14+HLA-DR+ cell vs CD14+HLA-DR− cell within 20 µm radius) was a prognostic factor independent of the densities of these cells, suggesting that proximity measurement of specific myeloid populations may improve the prognostic categorization compared with simpler, density-based analyses.

Each colorectal cancer is unique and arises through the complex interplay between endogenous and exogenous factors, neoplastic cells, and non-neoplastic cells.56 The current clinical practice recommendations for colorectal cancer are mainly based on disease stage,57 although there is an increasing need for more accurate prognostic and predictive parameters to guide personalized treatment. While our results support the prognostic significance of myeloid cells in colorectal cancer, further investigation is required to validate the findings and compare our assay to other methods for evaluating immune cell infiltrates, such as the T cell based Immunoscore that has been internationally validated using samples from more than 2600 patients.58

There are important limitations that need to be considered in the interpretation of these results. First, the study was based on TMAs. Myeloid immune infiltrates in a small tissue sample may not fully reflect overall tumor infiltration. However, multiple tumor cores were examined for most tumors (median 2), and we observed reasonably high core-to-core correlation, supporting the validity of the TMA approach. Second, while our multiplex immunofluorescence assay could detect more detailed myeloid cell phenotypes than standard immunohistochemistry, the number of markers was still limited. The best characterized human MDSC subsets include monocytic CD14+HLA-DR-CD15− cells and polymorphonuclear CD15+CD14− cells within the myeloid population, while less mature MDSCs are characteristically HLA-DR-CD33+.6 However, these phenotypic criteria do not enable unequivocal distinction of monocytic MDSCs from macrophages and monocytes, or polymorphonuclear MDSCs from mature neutrophils or eosinophils; functional suppression assays would be required to confirm the suppressive activity of the cells6 and are not compatible with archival FFPE material. Experimental studies evaluating the functional and biologic roles of various MDSC markers are also warranted, and orthogonal methods could be contrasted to protein expression based methods for MDSC-like cell detection. Third, data on cancer treatments were limited. Nonetheless, treatment decisions were likely based predominantly on disease stage, which was included in the multivariable Cox regression models. Future studies should investigate potential interactions between treatment and myeloid or monocytic cells in tumor tissues in relation to clinical outcomes.

The strengths of the study include the use of multiplex immunofluorescence combined with digital image analysis to phenotype granulocytic and monocytic myeloid cells in the tumor microenvironment. Our comprehensive tumor dataset enabled us to control for a large number of potential confounders, such as MSI status, CIMP status, KRAS, BRAF, and PIK3CA mutation status, and LINE-1 methylation. The study population was derived from cases in many hospitals across the USA, increasing the generalizability of the results compared with retrospective, single institution studies. We were able to control for selection bias due to tissue availability using the IPW method, as we could ascertain nearly all incident colorectal cancer cases in the cohort. Moreover, we used a stringent α level of 0.005 to reduce the possibility of false positive findings, thereby improving the likelihood that our findings are reproduced.34 However, additional validation in independent datasets will still be required to confirm the results.

In conclusion, myeloid cells exhibit strong, subset-specific prognostic significance in colorectal cancer, with mature CD14+HLA-DR+ and immature CD14+HLA-DR− monocytic phenotypes most notably showing opposite effects (favorable and unfavorable, respectively). These results support multimarker evaluation of myeloid immune infiltrates as a robust, quantitative prognostic tool in colorectal cancer.