Macrophages-based immune-related risk score model for relapse prediction in stage I-III non-small cell lung cancer assessed by multiplex immunofluorescence.

Background: Macrophages are critical players in regulating innate and adaptive immunity in the tumor microenvironment (TME). The prognostic value of macrophages and their heterogeneous phenotypes in non-small cell lung cancer (NSCLC) is still uncertain.
Methods: Surgically-resected samples of 681 NSCLC cases were stained by multiplex immunofluorescence to examine macrophage phenotypes as well as the expression levels of program death-ligand 1 (PD-L1) on them in both tumor nest and tumor stroma, including pan-macrophage (CD68+), M1 (CD68+CD163-), and M2 macrophages (CD68+CD163+). Various other immune cell markers, including CD4, CD8, CD20, CD38, CD66B, FOXP3, and CD133, were also evaluated. Machine learning algorithm by Random Forest (RF) model was utilized to screen the robust prognostic markers and construct the CD68-based immune-related risk score (IRRS) for predicting disease-free survival (DFS).
Results: The expression levels of CD68 were moderately correlated with the levels of PD-L1 (P<0.001), CD133 (P<0.001), and CD8 (P<0.001). Higher levels of CD68 (OR 1.03, 95% CI: 1.01-1.05, P<0.001) as well as M1 macrophage (OR 1.04, 95% CI: 1.01-1.06, P<0.001) indicated shorter DFS. Despite without statiscial significance, intratumoral M2 macrophage (OR 1.05, 95% CI: 0.99-1.10, P=0.081) was also associated with worse DFS. IRRS incorporating three intratumoral CD68-related markers and four intrastromal markers was constructed and validated to predict recurrence (high-risk group vs. low-risk group: OR 2.52, 95% CI: 1.89-3.38, P<0.001). The IRRS model showed good accuracy [area under the curve (AUC) =0.670, 0.709, 0.695, 0.718 for 1-, 3-, 5-year, and overall DFS survival, respectively] and the predictive performance was better than the single-marker model (area under the curve 0.718 vs. 0.500-0.654). A nomogram based on clinical characteristics and IRRS for relapse prediction was then established and exhibited better performance than the tumor-node-metastasis (TNM) classification and IRRS system (C-index 0.76 vs. 0.69 vs. 0.60, 0.74 vs. 0.67 vs. 0.60, 0.81 vs. 0.74 vs. 0.60 of the entire, training, testing cohort, respectively). Conclusions: Our study suggested close interactions between CD68 and other immune markers in TME, demonstrating the prognostic value of CD68 in relapse prediction in resectable NSCLC. 2022 Translational Lung Cancer Research. All rights reserved.

Introduction

Tumor progression relies on complex interactions between different components, such as adaptive and innate immune cells in the tumor microenvironment (TME), besides internal mechanisms in tumor cells (1). Consequently, evaluating the intrinsic and intricate structure in TME and identifying immune biomarkers with prognostic value contribute to elucidating the cellular and molecular mechanisms responsible for cancer development (2). Recently, the composition of immune contexture, including types [e.g., tumor-infiltrating lymphocytes (TILs)], quantity (e.g., cell density and proportion), and location (e.g., tumor islet or stroma), has been proven to play a pivotal role in the prognosis of non-small cell lung cancer (NSCLC) (3).

Given that prognosis varies considerably among NSCLC patients with similar pathological stage, evidently room for improvement still exists in the tumor-node-metastasis (TNM) classification (4). More recently, Galon et al. has proposed an impressive diagnostic improvement utilizing immune parameters (i.e., the density of CD3+ and CD8+ T cells), namely “Immunoscore”, as an essential supplement to the TNM stage model of colon cancer (5). In the field of NSCLC, through immune profiling of various biomarkers in TME, macrophages are of great interest.

Macrophages, mainly derived from bone marrow, are comprised of heterogeneous subsets. Specifically, tumor-associated macrophages (TAMs) are critical players in shaping the TME and capable to change their phenotypes based on ambient cytokine milieu (6). TAMs might have dual roles (i.e., anti-tumor and pro-tumor functions) because of two distinct subsets, the M1-polarized (classically activated) and M2-polarized (alternatively activated) populations (7). M1 macrophages, primarily induced by interferon-γ (IFN-γ), are considered to manifest anti-tumor effects by promoting Th1 responses (8). M2 macrophages, mainly activated by interleukin-4 (IL-4), can induce Th2 responses and show pro-tumor functions by immunosuppression and angiogenesis (9). Various markers have been implemented to identify different populations of TAMs, such as CD68 for pan-macrophages, human leucocyte antigen DR (HLA-DR) for M1, and CD163 for M2 macrophage, whereas no consensus concerning the optimal markers has been reached so far. For instance, CD68 is not exclusively expressed in TAMs, and other cell types like stromal cells might also express it in some way. Consequently, the amounts of macrophages might be overestimated, and the results can be biased. In addition to CD68 antibody specificity, the choice of different markers to phenotype M1 and M2 subsets may also contribute to their controversial prognostic roles reported by recent studies. In general, both positive (10), negative (11,12), or none (13,14) associations of TAMs with patients’ prognosis have been illustrated. The latest meta-analysis suggested that a higher density of TAMs in the tumor nest (TN) was associated with better overall survival (OS) while higher infiltrations of intrastromal TAMs predicted poor OS (15). In addition, they further found that higher infiltration levels of intratumoral M1 macrophages predicted better prognosis while a higher density of M2 macrophages in tumor stroma (TS) predicted worse prognosis.

In the current study, using the multiplex immunofluorescence (mIF) test for 681 stage I-III NSCLC primary tumor samples, we profiled the infiltration patterns of in situ macrophages according to broadly accepted markers of pan-macrophages (CD68), M1 macrophages (CD68+CD163-), and M2 macrophages (CD68+CD163+) as well as their functional status as assessed by the expression levels of programmed cell death-ligand 1 (PD-L1). Various other immune markers in TME, including CD4, CD8, CD20, CD38, CD66B, FOXP3, and CD133, were also evaluated. We further conducted the immune-related risk score (IRRS) model based on several CD68-related biomarkers identified as robust prognosticators by Random Forest (RF) algorithm. A nomogram system integrating the IRRS and clinicopathological features was also established as a quantitative tool to personalize risk assessment in predicting disease-free survival (DFS). We present the following article in accordance with the TRIPOD reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-21-916/rc).

Methods

Patient samples acquisition

Six hundred and eighty-one patients of stage IA to IIIB NSCLC who underwent radical resection (i.e., lobectomy/sub-lobectomy and lymph node dissection) were enrolled over the years 2009–2011 at the First Affiliated Hospital of Guangzhou Medical University, China. Additional inclusive criteria were: (I) all resected tissues and lymph nodes were confirmed by pathology finally and (II) sufficient resected tissues for mIF test. Patients were excluded if any of the following criteria were met: (I) multiple LC; (II) small cell lung cancer (SCLC) or non-invasive LC like lung adenocarcinoma (LUAD) in situ and minimally invasive LUAD; (III) diagnostic biopsy in pre-operation; and (IV) preoperative neoadjuvant therapy, as we previously described (16). The stage of primary tumors was classified based on the eighth edition of the TNM classification system (17). DFS is defined as the time from radical resection to local recurrence as measure by computed tomography (CT) scan. Patients received CT scan routinely every three months after the radical resection of lung cancer. Contrast-enhanced magnetic resonance imaging of head and radionuclide bone scanning were also done annually for detecting the recurrence (18). Ethics committee of the First Affiliated Hospital of Guangzhou Medical University approved this study (approval number: KLS-17-03), and all patients have signed the written informed consent for permitting mIF analyses of biological samples. The study was conducted following the Declaration of Helsinki (as revised in 2013) (19).

Multiplex immunofluorescence detection

The mIF staining was implemented at Genecast Biotechnology Co., Ltd. (Beijing, China). Cutting from formalin-fixed paraffin-embedded (FFPE) LC tissues, a four-µm thick slide was used for each panel test. The slides were dewaxed, rehydrated, and undergone epitope retrieval through boiling in Tris-EDTA buffer (pH=9; Klinipath #643901, the Netherlands) at 97 ℃ for 20 minutes (min). Subsequently, endogenous peroxidase was interdicted by incubation in Antibody Block/Diluent (PerkinElmer #72424205, USA) for 10 min, and protein was then blocked in 0.05% Tween solution for 30 min at 26 ℃. One antigen was tested in each round, including primary and secondary antibody incubation, tyramine signal amplification (TSA) visualization. Subsequently, the same procedure as before, i.e., targeting the next antibody after epitope retrieval and protein blocking. A total of 10 biomarkers, including CD66B, FOXP3, CD38, CD4, and CD20 in panel 1, and CD133, CD8, PD-L1, CD163, and CD68 in panel 2 were detected.

Primary antibodies for CD38 (ZM0422, clone SPC32, Zsbio, 1:400), CD66B (ab214175, polyclonal antibody, abcam, 1:50), CD8 (ZA-0508, clone SP16, Zsbio, 1:100), CD20 (ab9475, abcam, 1:50, Zsbio, 1:100), CD163 (ZM-0428, clone 10D6, PD-L1 (13684s, clone E1L3N, CST, 1:100), CD68 (ZM-0060, clone KP1, Zsbio, 1:100), FOXP3 (ab20034, clone 236A/E7, abcam, 1:100) were hatched for 1 h at 26 ℃, CD133 (ab19898, polyclonal antibody, abcam, 1:400) and CD4 (ZM0418, clone UMAB64, Zsbio, 1:200) were hatched for overnight at 4 ℃.

Anti-rabbit/mouse horseradish peroxidase (HRP) antibodies (Zsbio # PV-6002 or PV-8000) were utilized as the secondary antibody and hatched for 10 min at 37 ℃. TSA visualization was conducted by the opal seven-color mIF Kit (NEL797B001KT, PerkinElmer, USA), containing fluorophores (4',6-Diamidino-2-Phenylindole (DAPI)), Opal 570 (PD-L1 and CD4), Opal 520 (CD20 and CD163), Opal 650 (CD66B and CD133), Opal 620 (CD8), Opal 690 (CD68 and FOXP3), Opal 540 (CD38), and TSA Coumarin system (NEL703001KT, PerkinElmer, USA). Microwave treatment (MWT) was carried out to remove the TSA-antibody complex by Tris-EDTA buffer (pH=9; Klinipath #643901, the Netherlands) at 97 ℃ for 20 min, after targeting for all of the five antigens in each panel. All the slides were counterstained with DAPI for 5 min and were encapsulated in Antifade Mounting Medium (NobleRyder #I0052, China), prepared for imaging. Fresh whole-tissue slides cut from normal human tonsils were enrolled in each staining batch as positive control and evaluated the experiment’s reproducibility.

With the PerkinElmer Vectra (Vectra 3.0.5; PerkinElmer, USA), slides were scanned. Multispectral images were independent of spectral libraries built from single dyed tissue images for each antigen, employing the inform Advanced Image Analysis software (inForm 2.3.0; PerkinElmer, USA).

An algorithm is developed by training 10 to 15 typical multispectral images in batch analysis. Subsequently, tissue and cell segmentation were implemented by the algorithm. An experienced pathologist determined the appropriate positive threshold X for each biomarker. We defined X, 2X, 3X as the threshold of low fluorescence intensity (+), median fluorescence intensity (++), and high fluorescence intensity (+++), respectively. Histochemistry score (H-score) was calculated with the formula of H-score = (+++)%×3+ (++)%×2+ (+)%×1, as described by our previous published article (16).

Correlation analyses and prognostic analyses of CD68-related biomarkers

Using the Prism software (version 8.0), Spearman correlation coefficients were generated to investigate the associations between CD68-related biomarkers and the other biomarkers, including CD4, CD8, CD20, CD38, CD66B, FOXP3, and CD133. Heatmap and scatter plots were used to visualize the correlations. Since the data did not comply with normal distribution, the Kruskal-Wallis (K-W) test was utilized to compare the proportion of CD68-related biomarkers in TN and TS within the pathological slices of patients with or without relapse, and for comparing the proportion of CD68-related biomarkers in TN and TS among different clinical stage.

To study the independent prognostic value of the 108 immune biomarkers both in TN and TS, univariate and multivariate Cox analyses with covariates of age, sex, T stage, N stage, vascular cancer embolus, and number of lymph nodes resection, were performed using the SPSS software (version 25.0). The optimal cut-off values were determined based on the X-tile software (version 3.6) to divide the patients into high and low proportion groups (20). The prognostic significance of 26 CD68-related biomarkers was also evaluated by the Kaplan-Meier method and two-sided log-rank test.

Establishment of the CD68-based immune-related risk score by random forest algorithm

As a supervised learning approach based on random vectors of features, RF is a machine learning algorithm combining multiple tree predictors, which is capable of ranking the predictive ability of each variable and constructing a predictive model (21). We divided the entire cohort (n=681) into the training (n=477) and testing (n=204) cohort. Using the “randomForest” package (version 4.6-14) in R software (22), all the significant biomarkers identified by multivariate Cox regression (P<0.05) were inputted into the RF decision tree model. The optimal number of trees grown (ntree) was determined by the actual lowest mean squared error (MSE). Since a higher Mean Decrease Accuracy (%IncMSE) represents higher variable importance (23), biomarkers with a %IncMSE of above four were taken forward for further analysis. Ultimately, based on the selected biomarkers, the IRRS model was constructed by multiplication of the proportion and their regression coefficients originated from the multivariate Cox regression method:

where HR stands for the HR for biomarkers and proportion is their respective proportion. For investigating the prognostic significance of IRRS, we conducted multivariate Cox proportional hazards regression analysis adjusting for sex, age, T stage, N stage, number of lymph node resection, and vascular cancer embolus in the training cohort. According to the optimal cut-off value of the risk score, NSCLC patients in the training cohort were divided into high-IRRS and low-IRRS groups. Kaplan-Meier survival curves with log-rank test were exerted to contrast the DFS between the low-IRRS and high-IRRS groups. The receiver operating characteristic (ROC) curve and time-dependent ROC curve were then applied to assess the predictive accuracy of the IRRS model. Further verification of the robustness of IRRS was proceeded in the testing cohort and the entire cohort, utilizing the same formula and optimal cut-off value of the IRRS model in the training cohort.

Construction and validation of nomogram algorithm for relapse prediction

To provide the visualized risk prediction, we formulated a nomogram integrating IRRS and clinical characteristics based on multivariate Cox regression analysis results in the training cohort. Calibration plots were applied to estimate the predictive ability of the nomogram of 1-, 3-, and 5-year DFS. Nomograms and calibration plots were generated using the “rms” R package (version 6.2-0). In addition, with the “rmda” package (version 1.6) in R software, decision curve analysis (DCA) was employed for benefits evaluation and comparison of the nomogram model with standard TNM stage system and the IRRS model. Likewise, validation of the nomogram algorithm was applied in the testing cohort and the entire cohort.

Statistical analysis

Based on “survival” package (version 3.2-10) in R software, all survival curves were generated by the Kaplan-Meier method, and the two-sided log-rank test was used for disparity assessment of survival distributions. The chi-square test was applied respectively for rate comparison between two and multiple groups. Mann-Whitney U test and K-W test was exerted for the non-normal distribution data. Statistical analysis was performed using R (version 4.0.4, The R Foundation for Statistical Computing), SPSS (version 25.0, SPSS Inc), and Prism software (version 8.0, GraphPad Prism). All statistical tests were two-sided, and P value <0.05 was considered statically significant.

Results

Baseline characteristics of patients

A total of 681 eligible stage I to stage III NSCLC cases, with a median age of 60, were enrolled. Lung adenocarcinoma (479, 70.5%) and stage I to stage II patients (486, 74.8%) accounted for the majority of the cases. During the follow-up time, 282 patients relapsed (41.4%), of which 185 had vascular tumor embolus (VTE). Compared with patients without relapse, male patients (P=0.001) with more advanced TNM stage (P<0.001) have a higher tendency of recurrence. Additionally, patients with occurrence also tend to be with VTE (P<0.001) and visceral pleural invasion (VPI) (P=0.015) ().

Table 1

Clinicopathologic characteristics of the included patients

Characteristics	Without occurrence (n=399)	With occurrence (n=282)	Mann-Whitney U test (P value)
Age (%)
≤60, years	216 (54.1)	144 (51.1)	0.073
>60, years	183 (45.9)	138 (48.9)
Mean ± SD	61±11	58±11	–
Gender (%)			0.001
Male	212 (53.1)	186 (66.0)
Female	187 (46.9)	96 (34.0)
Histology (%)			0.178
Lung adenocarcinoma	289 (72.4)	190 (67.4)
Squamous cell lung cancer	84 (21.1)	66 (23.4)
Others	26 (6.5)	24 (8.5)
Pathology stage (%)			<0.001
IA	113 (28.3)	30 (10.6)
IB	152 (38.1)	50 (17.7)
IIA	49 (12.3)	56 (19.9)
IIB	15 (3.8)	21 (7.4)
IIIA	55 (13.8)	106 (37.6)
IIIB	1 (0.3)	2 (0.7)
T stage (%)			<0.001
T1	201 (50.4)	109 (38.7)
T2	134 (33.6)	80 (28.4)
T3	47 (11.8)	54 (19.1)
T4	17 (4.3)	35 (12.4)
N stage (%)			<0.001
N0	298 (74.7)	112 (39.7)
N1	51 (12.8)	99 (35.1)
N2	36 (9.0)	54 (19.1)
Vascular tumor emboli (%)			<0.001
Yes	230 (57.6)	185 (65.6)
No	169 (42.4)	97 (34.4)
Visceral pleural invasion (%)			0.015
PL0	171 (42.9)	91 (32.3)
PL1	200 (50.1)	172 (61.0)
PL2	28 (7)	19 (6.7)
Dissected lymph nodes (%)
0–14	100 (25.1)	83 (29.4)	0.633
≥15	285 (71.4)	182 (54.5)
Mean ± SD	20±10	21±10	–

CD68-related biomarker and its relationship with other biomarkers and clinical data

The proportion of immune biomarkers were proven to be non-normal distributing by the normality test. A total of 26 CD68-related biomarkers, including CD68 (POS), CD68+, CD68++, CD68+++, H-score (CD68), CD68+PDL1+, CD68+PDL1−, CD68+CD163+, CD68+CD163−, CD68+CD163+PDL1+, CD68+CD163+PDL1−, CD68+CD163-PDL1+, and CD68+CD163-PDL1−, and 92 other biomarkers both in TN and TS were tested and taken into analysis. Spearman rank correlation test () presented close relationships between the expression levels of intratumoral CD68+PDL1+ and intratumoral CD8 (r2=0.31, P<0.001) (), intratumoral CD68+CD163+ and PDL1 (r2=0.28, P<0.001) (), and intratumoral CD68 and intratumoral CD133 (r2=0.30, P<0.001) (). Highest infiltration levels of intratumoral CD68+++ macrophages, M2 macrophages, and M2 macrophages without expressing PDL1 were shown in stage IIB (). Significantly higher infiltration levels of pan-macrophages, M1 macrophages, and M1 macrophages without expressing PDL1 in TN than TS were observed in patients with relapse but not in patients without relapse ().

Figure 1

Spearman’s rank correlation between CD68-related biomarkers and other biomarkers in tumor nest and tumor stroma. (A) Spearman’s rank correlation matrix and corresponding P value between CD68-related biomarkers and various biomarkers in tumor and stroma. *, P<0.05; **, P<0.01; ***, P<0.001. (B) Merged immunofluorescence images show the expression profile of three biomarkers, including PDL1, CD8, and CD68. The dotted line graphs illustrate the correlations between (C) intratumoral CD68+PDL1+ and intratumoral CD8+CD133−; (D) intratumoral CD68+PDL1+ and intratumoral CD8 (POS). (E) Merged immunofluorescence images show the expression profile of three biomarkers, including CD68, CD163, and PDL1. The dotted line graphs illustrate the correlations between (F) intratumoral CD68+CD163+ and intratumoral PDL1 (POS); (G) intrastromal CD68+CD163+ and intrastromal PDL1 (POS). (H) Merged immunofluorescence images show the expression profile of two biomarkers, including CD68 and CD133. The dotted line graphs illustrate the correlations between (I) intratumoral CD68+CD163-PDL1− and intrastromal CD133 (POS); (J) intratumoral CD68+ and intratumoral CD133 (POS). Magnification of B,E,H: ×200.

Figure 2

Disparities of CD68-related biomarkers and their relationship with clinical data. Disparity of the proportion of CD68-related biomarkers in tumor nest (A) and tumor stroma (B) among different clinical stage. Disparity of the proportion of CD68-positive, CD68+CD163−, and CD68+CD163− PDL1-biomarkers in tumor nest and stroma in patients with (C-E) and without relapse (F-H). *, P<0.05; **, P<0.01. ns, not significant.

The prognostic values of CD68-related biomarkers

Univariate Cox regression analysis demonstrated that higher infiltration levels of intratumoral pan-macrophages and M1 macrophages were associated with significant shorter DFS [odds ratio (OR) 1.03, 95% confidence interval (CI): 1.01–1.04, P=0.004 for CD68 (POS); OR 1.04, 95% CI: 1.01–1.07, P=0.009 for CD68+; OR 1.03, 95% CI: 1.01–1.05, P=0.004 for CD68+CD163−]. Despite without statiscial significance, intratumoral M2 macrophage (OR 1.03, 95% CI: 0.99–1.10, P=0.081) was also associated with worse DFS. Moreover, these macrophages tended to be the subtypes without expressing PDL1 (OR 1.03, 95% CI: 1.01–1.04, P=0.003 for CD68+PDL1−; OR 1.03, 95% CI: 1.01–1.05, P=0.012 for CD68+CD163−PDL1−). The H-score of intratumoral CD68 also exhibited similar results (OR 1.01, 95% CI: 1.00–1.02, P=0.010). None of the intrastromal CD68-related biomarkers indicated significant relationships with DFS, implying the pivotal role of macrophages in TN (). We further cut the expression levels of CD68-related biomarkers into high- and low- value subtypes, and Kaplan-Meier curves were plotted to evaluate their prognostic significance (). Multivariate Cox regression analysis with adjustments for age, sex, T stage, N stage, VTE, and number of lymph nodes resection further supported that pan-macrophages (OR 1.03, 95% CI: 1.01–1.05, P<0.001) and M1 macrophages (OR 1.04, 95% CI: 1.01–1.06, P<0.001) rather than M2 macrophages (OR 1.05, 95% CI: 0.99–1.10, P=0.081) were independent prognosticators of worse DFS. Interestingly, we found that the prognostic impacts of CD68 varied from low fluorescence intensity (+) to high fluorescence intensity (+++). For instance, the prognosis effects of intratumoral CD68++ (OR 1.10, 95% CI: 1.02–1.19, P=0.015) were more significant than CD68+ (OR 1.05, 95% CI: 1.02–1.09, P=0.002) and CD68+++ (OR 1.07, 95% CI: 1.00–1.15, P=0.053), implying that the fluorescence intensity may represent the different developmental or functional states of macrophages (). Finally, a total of 29 biomarkers both in TN and TS, including 7 CD68-related biomarkers and 22 other biomarkers, which were independent prognosticators of DFS from multivariate Cox regression analysis, were selected as the candidate factors (Figure S1).

Figure 3

The prognostic values of CD68-related biomarkers. (A) Univariate Cox regression analysis and (B) multivariate Cox regression analysis evaluating the prognostic significance of CD68-related biomarkers for disease-free survival. Kaplan-Meier curves illustrate the prognostic implications of the proportion of CD68-related biomarkers (high vs. low) within tumor nest (C-M) and tumor stroma (N,O) for disease-free survival.

CD68-based immune-related risk score construction and its validation

Stratified sampling was utilized to divide 681 cases into a testing cohort (n=204, 20%) and a training cohort (n=477, 70%). The RF algorithm, with optimal parameter settings (mtry =3, ntree =368), was utilized to screen for the most important variables according to the expression levels of 29 candidate factors (). The importance of variable is ranked based on the %IncMSE, and the higher the %IncMSE, the more important a prognosticator is. Seven most significant predictors (i.e., %IncMSE >4), including intrastromal CD8+, intrastromal CD8+CD133−, intrastromal CD4+FOXP3+, intratumoral CD68+PDL1−, intrastromal CD4++, intratumoral CD68++, and intratumoral CD68+CD163− were identified and selected (). Subsequently, the IRRS of each patient was calculated using the formula: IRRS = ln (0.875) * (intrastromal CD8+)% + ln (0.945) * (intrastromal CD8+CD133−)% + ln (0.689) * (intrastromal CD4+FOXP3+)% + ln (1.035) * (intratumoral CD68+PDL1−)% + ln (0.836) * (intrastromal CD4++)% + ln (1.103) * (intratumoral CD68++)% + ln (1.035) * (intratumoral CD68+CD163−)%.

Figure 4

Construction the immune-related risk score (IRRS) by random forest (RF) model in the training cohort. (A) The importance of various variables within tumor nest and tumor stroma, evaluated by the value of Mean Decrease Accuracy (%lncMSE). (B) The importance feature map demonstrated the first rank of CD68-related biomarker in tumor nest. (C) Kaplan–Meier estimate of the disease-free survival based on IRRS (high IRRS group vs. low IRRS group) in the training cohort. (D,E) Time-dependent receiver operator characteristic (ROC) curves and area under curve (AUC) values indicated the accuracy of the IRRS model. (F) Comparison of the predictive performance between IRRS and single biomarkers using ROC curve analysis. (G) Immune landscape of two IRRS subgroups (high vs. low) illustrates distinct cellular characteristics. Infiltration disparities of two immune-related risk score subgroups (high vs. low) in tumor nest (H) and tumor stroma (I). Disparities of patients’ clinical characteristics of two IRRS subgroups (high vs. low) (J-P). IRRS, immune-related risk score; RF, random forest; %lncMSE, Mean Decrease Accuracy; ROC, receiver operator characteristic; AUC, receiver operator characteristic; ns, not significant; TPR, true positive rate; FPR, false positive rate.

With an optimal cut-off value (IRRS =−0.03), the Kaplan-Meier curve suggested patients with high IRRS were associated with a worse DFS than patients with low IRRS (median DFS survival: 899 vs. 2,415 days, P<0.001) in the training cohort (). Both univariate (OR 2.67, 95% CI: 2.03–3.60, P<0.001) and multivariate Cox regression analysis (OR 2.80, 95% CI: 1.87–3.98, P<0.001 in the training cohort; OR 2.40, 95% CI: 1.37–4.22, P=0.002 in the testing cohort; OR 2.52, 95% CI: 1.89–3.38, P<0.001 in the entire cohort) revealed an independently predictive role of IRRS in DFS of NSCLC (). The IRRS system also performed well in differentiating between patients with high-risk or low-risk recurrence (AUC =0.718, 95% CI: 0.668–0.768) (). We further implemented the time-dependent ROC curve analysis to evaluate its predictive performance, resulting in relatively high prediction accuracy (AUC =0.670, 0.709, 0.695 of 1-, 3-, and 5-year DFS survival, respectively) (). The AUC values exhibited that the prediction performance of the IRRS model was superior to a single marker (0.718 vs. 0.500–0.654) ().

Table 2

Construction and validation of immune-related risk score for predicting disease-free survival of non-small cell lung cancer

Variables	Univariate analysis			Multivariate analysis
				Training cohort			Testing cohort			Entire cohort
	HR (95% CI)	P value		HR (95% CI)	P value		HR (95% CI)	P value		HR (95% CI)	P value
Age	1.01 (1.00, 1.02)	0.055		1.02 (1.01, 1.04)	0.007		1.00 (0.98, 1.03)	0.908		1.02 (1.00, 1.03)	0.013
Sex
Male	Ref.			Ref.			Ref.			Ref.
Female	0.70 (0.54, 0.90)	0.005		0.73 (0.52, 1.04)	0.080		0.79 (0.44, 1.41)	0.419		0.76 (0.56, 1.02)	0.066
T stage
T1	Ref.			Ref.			Ref.			Ref.
T2	1.05 (0.77, 1.41)	0.770		0.75 (0.50, 1.14)	0.174		0.90 (0.45, 1.81)	0.766		0.80 (0.57, 1.13)	0.206
T3	1.64 (1.16, 2.32)	0.005		1.51 (0.93, 2.46)	0.098		1.52 (0.68, 3.40)	0.310		1.45 (0.97, 2.18)	0.069
T4	2.85 (1.93, 4.22)	<0.001		1.64 (0.94, 2.84)	0.081		2.87 (1.14, 7.22)	0.025		1.88 (1.18, 2.99)	0.007
N stage
N0	Ref.			Ref.			Ref.			Ref.
N1	2.85 (2.03, 4.00)	<0.001		2.38 (1.44, 3.95)	0.001		2.57 (1.21, 5.46)	0.014		2.37 (1.57, 3.58)	<0.001
N2	3.40 (2.56, 4.52)	<0.001		2.87 (1.92, 4.30)	<0.001		5.91 (2.99, 11.68)	<0.001		3.43 (2.45, 4.81)	<0.001
Visceral pleural invasion
PL0	Ref.			Ref.			Ref.			Ref.
PL1	1.45 (1.11, 1.90)	0.006		1.16 (0.78, 1.71)	0.471		1.01 (0.54, 1.87)	0.981		1.12 (0.81, 1.53)	0.498
PL2	1.26 (0.76, 2.11)	0.370		1.27 (0.64, 2.52)	0.496		1.26 (0.29, 5.59)	0.759		1.24 (0.68, 2.28)	0.487
Vascular tumor emboli
No	Ref.			Ref.			Ref.			Ref.
Yes	2.05 (1.58, 2.65)	<0.001		1.45 (0.99, 2.12)	0.054		1.56 (0.83, 2.92)	0.170		1.47 (1.06, 2.02)	0.019
Resected lymph nodes	1.00 (0.98,1.01)	0.613		0.98 (0.96, 1.00)	0.017		0.98 (0.95, 1.01)	0.123		0.98 (0.97, 0.99)	0.005
Immune-related risk score
Low	Ref.			Ref.			Ref.			Ref.
High	2.67 (2.03, 3.60)	<0.001		2.80 (1.87, 3.98)	<0.001		2.40 (1.37, 4.22)	0.002		2.52 (1.89, 3.38)	<0.001

HR, hazard ratio.

The same formula and cutoff value of IRRS model were subsequently applied to the testing cohort and the entire cohort. Consistent with the findings from the training cohort, the high-IRRS subgroup had a significantly lower DFS compared with the low-IRRS subgroup (median DFS survival: 1,073 vs. 2,355 days, P<0.001 in the testing cohort; median DFS survival: 1,022 vs. 2,415 days, P<0.001 in the entire cohort). ROC curve further confirmed its relatively great sensitivity and specificity in DFS prediction (AUC =0.676, 95% CI: 0.595–0.757 in the testing cohort; AUC =0.704, 95% CI: 0.661–0.747 in the entire cohort) (Figures S2,S3).

Cellular and clinical characteristics of immune-related risk score

The discrepancies of cellular characteristics between the high-IRRS and low-IRRS groups in the training cohort were illustrated by a sweeping landscape (). K-W tests revealed the expression levels of intratumoral CD68-related biomarkers, including CD68+, CD68+PDL1−, CD68+CD163−, and CD68+CD163−PDL1−, were significantly higher in the high-IRRS group than the low-IRRS group (). However, the expression levels of intrastromal biomarkers, including CD4+, CD4++, CD4+++, CD38+, CD4+FOXP3+, CD4+FOXP3−, CD8+, CD8++, CD8+CD133−, and PDL1+ were significantly higher in the low-IRRS group than the high-IRRS group, implying immune-activated microenvironment in the low-IRRS subtype (). Chi-square tests showed that the clinical stage (P<0.05) and T stage (P<0.05) were more advanced in the high-IRRS group compared with the low-IRRS group whereas the N stage, VPI, and VTE were not ().

In the testing dataset, analogously, the expression levels of CD68++, CD68+PDL1−, CD68+CD163−, and CD68+CD163-PDL1− were significantly higher while CD4+, CD4++, CD4+++, CD38+, CD4+FOXP3+, CD4+FOXP3−, CD8+, CD8++, and CD8+CD133− were significantly lower in the high-IRRS group than the low-IRRS group (Figure S2C,S2D). The T stage (P<0.05) was more advanced, and patients’ age (P<0.05) was older in the high-IRRS group than the low-IRRS group (Figure S2E-S2K).

Establishment of relapse prediction system for individual NSCLC patients

To construct a quantitative tool and personalize risk assessment to predict DFS, a nomogram integrating IRRS and clinicopathological characteristics, including sex, age, VTE, TNM stage, and the number of resected lymph nodes, was proposed utilizing the training dataset. The IRRS contributed slightly less to the risk points than the clinical stage, while similar to the prognostic effects of the number of resected lymph nodes (). The calibration plots fitted well to the ideal curve in 1-year DFS, whereas it was not observed in 3- or 5-year DFS (). The C-index of nomogram system (C-index =0.74, 95% CI: 0.70–0.78) was higher than the standard TNM stage (C-index =0.67, 95% CI: 0.64–0.71) and the IRRS model (C-index =0.60, 95% CI: 0.57–0.63), suggesting the nomogram algorithm improved the effectiveness of predictive model (). We further implemented the DCA, indicating that the nomogram added more benefit to predict patients’ DFS than IRRS system and TNM stage (). Great calibrations were also illustrated for the probability of recurrence in 1-year DFS in the testing and entire cohort. Analogously, the prognostic nomogram showed better predictive capability than the standard TNM stage and the IRRS model (C-index =0.81 vs. 0.74 vs. 0.60 in the testing cohort; C-index =0.76 vs. 0.69 vs. 0.60 in the entire cohort) (, Figures S4,S5). Based on the above-mentioned findings, the IRRS model might be employed to reinforce the prognostic ability of the TNM approach.

Figure 5

Construction of nomogram for predicting the disease-free survival in the training cohort. (A) The nomogram was constructed with the immune-related risk score in the training cohort. (B) Forest plot of the multivariate Cox regression analysis. Calibration plot and decision curve analysis of the nomogram in terms of agreement between the predicted and observed (C,F) 1-, (D,G) 3- and (E,H) 5-year outcomes.VCE, vascular cancer embolus; IRRS, immune-related risk score; LNs, lymph nodes; VPI, visceral pleural invasion; AIC, Akaike information criterion; DFS, disease-free survival.

Table 3

Harrell’s concordance indexes of TNM stage, immune-related risk score, and nomogram system

Cohort	TNM stage (95% CI)	IRRS (95% CI)	Nomogram (95% CI)
Training	0.67 (0.64, 0.71)	0.60 (0.57, 0.63)	0.74 (0.70, 0.78)
Testing	0.74 (0.68, 0.80)	0.60 (0.56, 0.64)	0.81 (0.76, 0.86)
Entire	0.69 (0.66, 0.73)	0.60 (0.54, 0.66)	0.76 (0.73, 0.79)

IRRS, immune-related risk score; CI, confidence interval.

Discussion

The high heterogeneity of TAMs, characterized as the variability of their constituents depending on clinical stage, cancer types, and intratumor heterogeneity, has been studied by numerous researchers, revealing their vital role in cancer development, progression, and recurrence (24-27). CD68, considered as the significant biomarker for quantifying macrophage amounts, is generally associated with prognosis and correlated with resistance to cancer treatment in solid tumors, especially lung cancer, which is the deadliest cancer worldwide (28). Based on the mIF approach, for the first time, we delved into the clinical significance of CD68-related biomarkers in TME of stage I to III NSCLC and appraised their correlations with other biomarkers. Further, an effective CD68-based IRRS system based on seven immune biomarkers was constructed by the RF algorithm. Moreover, integrating the IRRS, our nomogram model was able to provide a personalized prognosis tool for NSCLC patients.

Based on univariate Cox regression and survival analysis, the proportion of intratumoral and stromal CD68-related biomarkers were principally the risk factors for DFS in lung cancer, consistent with previous findings (29,30). Furthermore, the multivariate Cox regression illustrated the intratumoral proportion of CD68 (including CD68+, CD68++, CD68+PDL1−, CD68+CD163−, and CD68+CD163−PDL1−), alongside its H-score, were the independent prognostic indicators for poorer DFS. Characterized as one of the main contributors of tumorigenesis and a major component of the malignant microenvironment, CD68+ macrophages are involved in a wide variety of mechanisms related to tumorigenesis (31). During the early formation of NSCLC, aggregation of TAMs took place nearby the tumor cells, thereby promoting the invasiveness and epithelial-mesenchymal transition (EMT). Likewise, they induce a potent regulatory T cell response that restrains adaptive antitumor immunity (32). This correlation is well demonstrated in our study, that M2 macrophages were significantly correlated with PD-L1 positive cells in TS (r2=0.31, P<0.001) and TN (r2=0.28, P<0.001) respectively, thereby directly inhibiting the antitumor response through cell-cell interaction with cytotoxic T cells and the production of immunosuppressive cytokines such as transforming growth factor-β (TGF-β) and IL-10 (33).

Besides, we found a close relationship between intratumoral CD68 and intratumoral CD133 (r2=0.30, P<0.001), supporting the bidirectional crosstalk between cancer stem cells and TAMs, which contributes to a pro-tumorigenic and immunosuppressive contexture (34), leading to the poorer prognosis of the patients with higher pan-macrophages infiltration in TME. In fact, our findings supported that the proportion of pan-macrophages, mainly M1 macrophages, were higher In TN among the recurrence group, while this disparity was not found in the non-recurrence group, further supporting the vital role for CD68 in the progression and recurrence of lung cancer. Interestingly, PD-L1-positive pan-macrophages, especially M1 macrophages, were prone to associate with longer DFS according to the multivariate Cox regression analysis, which is similar to the findings of a previous study (35). Simultaneously, the level of intratumoral PDL1-positive pan-macrophages was found to associate with the level of CD8-positive T cells significantly (r2=0.31, P<0.001) in our study, suggesting a connection between high PD-L1 and immune-inflamed tumors, namely “hot tumors.” Yet, the innate mechanisms remained uncertain and warranted further investigations.

However, not all the studies were consistent. An improved survival in NSCLC patients was confirmed in a multiplex immunohistochemistry (mIHC) based research by Rakaee et al. aimed at investigating the prognostic value of macrophages (36). Specifically, this study suggested an improved NSCLC-specific survival with high M1 and M2 infiltration levels in both TN and TS (36). The source of the discrepancies may originate from three aspects. Firstly, due to the methodological variation between both studies, inconsistent antibodies were used for targeting biomarkers, leading to the possible disparity of precise definition of the immunophenotyping of macrophages (37). Secondly, given the assay of different methodologies, namely mIHC and mIF, mIHC has a limited dynamic range of most chromogenic substrates, and only very few can be combined for studying marker co-expression. In contrast, a more extensive linear dynamic range can be detected using mIF, which facilitates the comprehensiveness and precision when quantifying various biomarkers on a single slide (37). Thirdly, it was worth mentioning that the survival outcomes of the study by Rakaee et al. differed from ours, in which disease-specific survival (DSS) was used as the endpoint with longer follow-up time (median follow-up was 86 months) (36). In this case, we hypothesize that although more macrophages infiltration in TME was associated with tendency of cancer recurrence at an earlier time, they might have a positive impact on the survival of patients who have relapsed, as evidence suggested that TAMs could improve the response to chemotherapy of LC patients and extent their survival (24,38,39). Albeit chemotherapy was traditionally considered immunosuppressive due to potential immune cells depletion, chemotherapy could promote immunogenic cell death in cancer and conduce type-1 IFN signaling, thereby enhancing the recruitment of antigen presenting cells to the TME and facilitating efficient T cell priming (40,41).In the field of epidermal growth factor receptor (EGFR)-mutant LC patients, the prognostic value of macrophages has been a research focus in recent years. Saxon et al. reported that epithelial NF-κB signaling sustained carcinogenesis in EGFR-mutant LC by recruiting pro-tumorigenic macrophages (42). Wang et al. further found that EGFR oncogene–dependent progression was correlated with the expansion of alveolar macrophages and the presence of activated signature of macrophages indicated unfavorable OS in patients receiving resection for EGFR-mutant LUAD (43). Yin and colleague showed that nanomedicine remodeling the tumor microenvironment (e.g., M2-macrophage repolarization) could reverse the resistance of EGFR-tyrosine kinase inhibitor and enhance treatment outcomes in mouse model (44). As for the effect of how neoadjuvant therapy influences TAMs within the tumor environment, Parra et al. (45) discovered that the density of CD68+ TAMs was higher in neoadjuvant chemotherapy (NCT) NSCLC than in non-NCT-NSCLC cases. Besides, a recent study conducted by Gaudreau et al. (46) also revealed that patients with NSCLC who received NCT followed by surgery, as compared with patients who received upfront surgery (US), had overall higher levels of immune infiltration, including higher densities of CD68+TAMs in their tumors. As suggested by Blankenstein (47), in the tumor epithelial and stromal compartments, the activation of TAM class M1 and helper T cells (CD3+ CD4+) mediate tumor suppression factors in NCT-treated NSCLC patients. Considering the highly functional heterogeneity of macrophages in TME, extensive efforts are still required to corroborate our hypothesis.

Meanwhile, several limitations in previous reports should also be noted, such as small sample size and homogeneous cohorts (i.e., utilizing a specific TNM stage), which may lead to the contradictory findings. Moreover, the wide variation of methods used in evaluating the infiltrating patterns of TAMs in previous studies made it difficult to draw a definitive conclusion concerning their associations to prognosis. Despite that immunohistochemistry (IHC) is a helpful tool in diagnostic settings and has been widely used for decades, some drawbacks limit its ability to assess the complex components in TME. For instance, the subjectivity in the interpretation of IHC stain by pathologists may significantly affect the reproducibility of CD68 expression levels (48). Additionally, although multiplex IHC (mIHC) can simultaneously identify colocalized markers, it is limited by quantitative assessment (i.e., appraise expression levels as positive vs. negative simply or a semiquantitative H-score) (49). The mIF approach is an emerging and powerful approach, based on the principle that various protein targets can be dyed apiece by a particular antibody labeled with a disparate fluorophore (37). Given that the fluorophores have a sizeable dynamic scope and be captured in situ by the multispectral microscope, the mIF tool can characterize cells phenotypically and facilitate quantitative and spatial analysis.

As evidence from current and previous studies suggests that CD68-associated biomarkers have significant implications as prognostic markers for clinical outcomes of NSCLC, in order to access the value for cancer relapse and the underlying therapeutic implications, it is essential to develop a helpful tool for prognosis prediction. In this study, the RF algorithm was implemented, and the proportion of intratumoral CD68+PDL1−, CD68++, CD68+CD163−, and intrastromal CD8+, CD8+CD133−, CD4+FOXP3+, CD4++ were selected and combined to construct an IRRS model for predicting DFS. This CD68-based prognostic model has not been reported and may represent a novel prognostic factor for NSCLC. The AUCs for the total and 1-, 3-, and 5-year DFS rates for the IRRS were 0.718, 0.670, 0.709, and 0.695 respectively in the training cohort (n=477). For providing a personalized scoring system for the prognosis of each NSCLC patient, we established a nomogram combining the IRRS and clinicopathological characteristics, including sex, age, vascular tumor emboli (VTE), TNM stage, and the number of resected lymph nodes). This nomogram system could further improve the stratification ability of the prognosis, as the C-index of which was 10.1% higher than the standard TNM stage (C-index =0.76 vs. 0.69) for the entire cohort (n=681).

Since our IRRS model presents prominent prognostic implications, we further probed into the characteristics of the immune landscape between the high and the low IRRS groups. The infiltration disparities of immune biomarkers could provide clues for the difference in recurrence and survival prognosis. The proportion of CD38+, CD4-related biomarkers (including CD4+, CD4++, CD4+++, and CD4), CD8-related biomarkers (including CD8+, CD8++, and CD8+CD133−) in TS were higher in the low IRRS group. Instead, the high IRRS group demonstrated a significantly higher proportion of CD68-related biomarkers infiltration (including CD68+, CD68++, CD68+PDL1−, and CD68+CD163−). Besides, high IRRS group was associated with a more advanced T stage and clinical stage, as a higher amount of CD68+ macrophages were positively correlated with tumor growth and metastasis, consistent with the conclusions of previous studies (24,50). These findings are reasonably consistent with the clinical practice. It is well studied that CD4+ and CD8+ T cells were inhibited by the maintenance of strong immunosuppressive TME due to T cell immune checkpoint ligands expressed or cytokines secreted by TAMs (51). In the context of the success of ADAURA (52) and IMPower 010 (53) that adjuvant therapies brought pronounced benefit for resected early-stage NSCLC cases, targeting TAMs is emerging as an attractive and effective strategy (ClinicalTrials.gov identifier: NCT05053750, NCT01765790, NCT00317603, etc.) for therapeutic intervention (54). During the course of therapy, the efficacy of macrophages depletion is constantly relying upon the function and enhanced recruitment of cytotoxic CD8+ T cells (55), allowing an effective antitumoral immune response. One of the most studied therapeutic targets, the macrophage colony stimulating factor-1 (CSF-1), was revealed to recruit inflammatory monocytes into the cancer site. Due to the presence of cytokines such as IL-10 and IL-13, they subsequently differentiate into TAMs, which contribute to the formation of a complex TME that promotes malignant transformation and cancer progression. After the inhibition of CSF-1, the infiltration of TAMs was abrogated, and the recruitment of CD8+ T cells was enhanced, thereby restricting tumor growth (56).

Several limitations need to be mentioned in our study. First, it is worth noting that our study was based on the percentage of immune biomarkers in pathological slices. Therefore, spatial variables, for example, the proximity effects between immune cells were not fully considered in our study. Second, due to the research design, we did not realize the validation of the above biomarkers and risk model on NSCLC patients’ DFS from external cohorts. Third, only cytological analyses were provided in the current study. Hence, multi-omics studies are needed to clarify the regulatory mechanism between macrophages and other components in the TME of NSCLC.

Nevertheless, considering the therapeutic and prognostic value of CD68, based on the mIF assay, the IRRS model and nomogram system were valuable signatures for their correlations with immune landscapes of NSCLC. Besides, the association between DFS and the CD68-based IRRS model in the comprehensive set of 681 patients indicated its powerful prognostic marker for NSCLC. Our study may help define the prognosis of NSCLC patients further, provide enlightenment for further exploration of the role of macrophages in TME, and effectively support the exploration of TAMs-related therapeutic approaches in NSCLC.

The article’s supplementary files as