Article Text

Original research
Immune suppressive microenvironment in liver metastases contributes to organ-specific response of immunotherapy in advanced non-small cell lung cancer
  1. Jia-Yi Deng1,2,
  2. Qing Gou3,
  3. Lingling Yang4,
  4. Zhi-Hong Chen1,
  5. Ming-Yi Yang1,2,
  6. Xiao-Rong Yang1,2,
  7. Hong-Hong Yan1,
  8. Xue-Wu Wei1,2,
  9. Jia-Qi Liu1,
  10. Jian Su1,
  11. Wen-Zhao Zhong1,
  12. Chong-Rui Xu1,
  13. Yi-Long Wu1 and
  14. Qing Zhou1
  1. 1Guangdong Lung Cancer Institute, Guangdong Provincial People's Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou, China
  2. 2School of Medicine, South China University of Technology, Guangzhou, China
  3. 3Department of Interventional Radiology, Guangdong Provincial People's Hospital (Guangdong Academy of Medical Sciences), Southern Medical University, Guangzhou, China
  4. 4Geneseeq Research Institute, Nanjing Geneseeq Technology Inc, Nanjing, China
  1. Correspondence to Professor Qing Zhou; gzzhouqing{at}126.com

Abstract

Background The liver is a frequent site of metastases and liver metastases (LM) correlate with diminished immunotherapy efficacy in non-small cell lung cancer (NSCLC). This study aimed to analyze whether tumor response to immunotherapy differs between pulmonary lesions (PL) and LM in NSCLC and to explore potential mechanisms through multiomics analysis.

Methods This observational longitudinal clinical cohort study included patients with NSCLC with LM receiving immunotherapy was conducted to evaluate organ-specific tumor response of PL and LM. We collected paired PL and LM tumor samples to analyze the organ-specific difference using whole-exome sequencing, RNA sequencing, and multiplex immunohistochemistry.

Results A total of 52 patients with NSCLC with LM were enrolled to evaluate the organ-specific response of immunotherapy. The objective response rate (21.1% vs 32.7%) and disease control rate of LM were lower than that of PL (67.3% vs 86.5%). One-third of patients showed mixed response, among whom 88.2% (15/17) presented with LM increasing, but PL decreasing, while the others had the opposite pattern (p=0.002). In another independent cohort, 27 pairs of matched PL and LM tumor samples from the same individuals, including six simultaneously collected pairs, were included in the translational part. Genomic landscapes profiling revealed similar somatic mutations, tumor mutational burden, and neoantigen number between PL and LM. Bulk-RNA sequencing showed immune activation-related genes including CD8A, LCK, and ICOS were downregulated in LM. The antigen processing and presentation, natural killer (NK) cell-mediated cytotoxicity and T-cell receptor signaling pathway were enriched in PL compared with LM. Multiplex immunohistochemistry detected significantly lower fractions of CD8+ cells (p=0.036) and CD56dim+ cells (p=0.016) in LM compared with PL. Single-cell RNA sequencing also characterized lower effector CD8+ T cells activation and NK cells cytotoxicity in LM.

Conclusions Compared with PL, LM presents an inferior organ-specific tumor response to immunotherapy. PL and LM showed limited heterogeneity in the genomic landscape, while the LM tumor microenvironment displayed lower levels of immune activation and infiltration than PL, which might contribute to developing precise immunotherapy strategies for patients with NSCLC with LM.

  • Tumor Microenvironment
  • Immunotherapy
  • Non-Small Cell Lung Cancer
  • Lung Neoplasms
  • Liver Neoplasms

Data availability statement

Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information.

http://creativecommons.org/licenses/by-nc/4.0/

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.

Statistics from Altmetric.com

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • Previous studies reported organ-specificity in the tumor microenvironment (TME) at different anatomical sites. However, the comparison of the TME between strictly paired pulmonary lesions (PL) and liver metastases (LM) in primary non-small-cell lung cancer (NSCLC) has not been explored due to the low accessibility of paired samples.

WHAT THIS STUDY ADDS

  • We concluded that LM presents an inferior organ-specific tumor response to immunotherapy compared with PL. This is the first study that used strictly paired NSCLC samples to explore organ-specific TME between PL and LM, from which we demonstrated that LM TME generally exhibits lower levels of immune infiltration and immune activation.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • These findings provide insights into the mechanisms underlying the organ-specific tumor response to immunotherapy and information for the development of precise immunotherapy strategies for patients with NSCLC with LM.

Background

Significant progress has been made in systemic therapy for advanced non-small cell lung cancer (NSCLC) in recent years. Several phase III trials have demonstrated the significant survival benefits of immunotherapy targeting immune checkpoints, including programmed death-1 (PD-1)/programmed death-ligand 1 (PD-L1) interaction, in patients with advanced NSCLC.1 The liver, which is richly supplied with portal venous and arterial blood supply, is one of the most frequent sites of tumor metastasis. It is also considered an immune-tolerant organ, characterized by T-cell anergy and immunosuppressive signals.2 Recent studies identified that liver metastases (LM) could cultivate an immune desert, which may be associated with inferior response to immunotherapy, as well as shorter overall survival (OS) and progression-free survival (PFS) in patients with advanced solid tumors.3 4 Noteworthily, a previous study revealed that LM could actively induce CD8+ T cells to undergo apoptosis following their interaction with macrophages, concluding that LM co-opt host peripheral tolerance mechanisms to cause acquired immunotherapy resistance through CD8+ T-cell deletion.3

Immune checkpoint inhibitors (ICIs) induce antitumor effects by reactivating exhausted T cells and thus restoring systemic antitumor immunity, in which the tumor microenvironment (TME) plays a vital role. The dominant forces that influence the rate or pattern of tumor growth include the host’s innate and adaptive immune response, tumor-cell-intrinsic properties, and the anatomical organ in which the tumor resides.5 Heterogeneous TME of various organs, called organ-specific TME, may potentially influence tumor growth and responses to anticancer treatment.6–8 The success of ICIs against both primary and metastatic aerodigestive malignancies may be a prime example of the importance of the site of the TME.9 Recently several studies reported this phenomenon where patients with metastatic urothelial carcinoma or hepatocellular carcinoma experienced mixed responses in different metastasis sites while receiving ICIs. This phenomenon in which tumors in one organ enlarge while those in another organ ameliorate or remain stable was called the organ-specific mixed response. These support that tumors in the liver, whether primary or metastatic, are more likely to resist immunotherapy as the liver is viewed as an immune-tolerant organ.6 10 These differences have generally been attributed to the genetic heterogeneity of the tumor cell itself, as well as organ-specific features of the TME. However, previous studies about organ-specific TME have always used animal models or non-strictly paired NSCLC tumor samples, which could not exclude interpatient heterogeneity.11 12

Our prior studies found that metastases in different anatomical locations may be associated with different clinical outcomes and local tumor responses to ICIs in NSCLC.13 Herein, we designed this study, which was based on strictly paired tumor samples from the same individuals, to further analyze whether tumor response to ICIs differs between pulmonary lesions (PL) and LM in NSCLC and to explore potential cellular and molecular mechanisms underlying this clinical phenomenon by multiomics analysis.

Methods

Clinical cohort design and participants

This observational longitudinal cohort study included patients with consecutive advanced NSCLC with LM receiving immunotherapy who were treated at Guangdong Provincial People’s Hospital between May 2017 and April 2021. The participants were enrolled from a prospectively managed database that collected information on their clinical characteristics, treatment, and clinical outcomes. The inclusion criteria were the following: (1) age ≥18 years; (2) diagnosis of NSCLC with LM; (3) receiving at least one course of ICIs; and (4) with at least one evaluable target lesions in both liver and lung at baseline and undergoing imaging for tumor response assessments every two to three courses of the treatment until the discontinuation of ICIs treatment. Patients with incomplete or missing data were excluded from the final analysis. The study investigators played no role in the design or implementation of the treatment protocols. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

An independent oncologist conducted organ-specific tumor response evaluation of clinical imaging in tumors located at the lung and liver, respectively. The organ-specific response criteria, according to the response evaluation criteria in solid tumors (RECIST) V.1.1 criteria, were applied. All lung lesions and liver metastatic lesions that met the RECIST V.1.1 criteria were included in the analysis. The organ-specific tumor response, including complete responses (CR; complete disappearance), partial response (PR; ≥30% reduction, taking as reference the baseline sum of diameters), progressive diseases (PD; ≥20% increase, taking as reference the smallest sum of diameters), or stable diseases (SD; neither CR, PR, nor PD), were determined for the lung and liver systems, respectively. Organ-specific response rate (OSRR) was defined as the percentage of patients who had CR or PR as the best response of the target lesions in each organ. Organ-specific disease control rate (OSDCR) was defined as the ratio of the number of patients reaching CR, PR, or SD as the best organ-specific tumor response to the total number of patients. PFS and OS were calculated from the time of therapy initiation.

Paired tumor samples collection

Strictly paired PL and LM tumor samples were obtained from the biospecimen bank of the Guangdong Lung Cancer Institute. All samples were collected within 30 min after resection and dissected into a 0.5 cm3 block. Samples were freshly frozen in liquid nitrogen, stored at –80°C, and/or paraffin‐embedded after decalcification and stored at room temperature. H&E-stained sections were assessed by a pathologist to confirm NSCLC content. Specimens that did not meet minimal tissue requirements for assessment were excluded. A total of 52 archived frozen or formalin-fixed, paraffin-embedded (FFPE) tissues from 26 patients were used. Another pair of PL and LM fresh tumor samples from Pa27 were prospectively collected. According to the timing of sampling, paired samples were divided into four groups: (1) PL−LM: PL samples were collected at radical operation, while LM samples were collected at the time of first postoperative recurrence; (2) PL=LM: PL and LM samples were simultaneously collected; (3) PL>LM: PL samples were collected anterior to LM; and (4) PL<LM: PL samples were collected posterior to LM. Informed patient consent was obtained through the biospecimen bank of the Guangdong Lung Cancer Institute for specimen use.

Whole-exome sequencing

Fresh frozen tissues from tumor samples were used for genomic DNA extraction with QIAamp DNA FFPE Tissue Kit (QIAGEN) following the manufacturer’s instructions. White blood cell DNA was sequenced together with tumor DNA samples to identify germline mutations. The DNA quality was assessed by NanoDrop 2000 (Thermo Fisher Scientific) and the quantity was measured by dsDNA HS Assay Kit (Life Technologies) on Qubit V.2.0. Extracted tumor genomic DNA was fragmented into 300~350 bp using the Covaris M220 instrument (Covaris). Sequencing libraries were prepared with KAPA HyperPrep kit (KAPA Biosystems) with optimized protocols. DNA libraries from different samples were marked with unique indices during library preparation and up to 2 µg of different libraries were pooled together for targeted enrichment. The enriched libraries were sequenced on HiSeq 4000 NGS platforms (Illumina) to coverage depths of 200×. Trimmomatic was used for FASTQ file quality control. GATK V.3.4.0 was applied to detect germline mutations from blood control samples. VarScan V.2 was employed for detection of somatic mutations. Annotation was performed using ANNOVAR using the hg19 reference genome. TMB was defined as the total single nucleotide variant (SNV) counts in coding regions14 and neoantigen was analyzed by NeoPredPipe (Neoantigen Prediction Pipeline).15

Bulk-RNA sequencing and analysis

Total RNA was extracted from samples using the RNeasy FFPE kit (QIAGEN). Ribosomal RNA was depleted using RNase H followed by library preparation using the KAPA Stranded RNA-seq Kit with RiboErase (HMR) (KAPA Biosystems). Library concentration was determined by the KAPA Library Quantification Kit (KAPA Biosystems), and library quality was accessed by the Agilent High Sensitivity DNA kit on a Bioanalyzer 2100 (Agilent Technologies). The library was then sequenced on an Illumina HiSeq NGS platform (Illumina). The CIBERSORT16 algorithm was employed to quantify the proportions and distributions of tumor-infiltrating immune cells based on the RNA-sequencing data. Gene set enrichment analysis (GSEA) was performed using GSEA software (V.4.1.0).17

Multiplex immunohistochemistry

Multiplex immunohistochemistry (mIHC) staining and multispectral imaging were obtained to identify the immune cell distribution in the TME using PANO 7-plex IHC kit (Panovue, Beijing, China). Staining for CD56 (CST3576, CST, Massachusetts, USA), CD68 (CST76437, CST, Massachusetts, USA), HLA-DR (ab92511, Abcam, Massachusetts, USA), CD8A (CST70306, CST, Massachusetts, USA) and panCK (CST4545S, CST, Massachusetts, USA) was performed sequencing using primary antibodies. Different primary antibodies were sequentially applied, followed by horseradish peroxidase-conjugated secondary antibody incubation and tyramide signal amplification. The slides were microwave heat-treated after each tyramide signal amplification (TSA) operation. Nuclei were stained with 4′−6′-diamidino-2-phenylindole (DAPI, Sigma-Aldrich) after all the human antigens had been labeled. These markers were used to simultaneously detect T cells (CD8A), natural killer (NK) cells (CD56), macrophages (CD68/human leukocyte antigen-DR (HLA-DR)), and tumor cells (pan cytokeratin (panCK)). All tumor samples were divided into the tumor parenchyma area (tumor cell-enriched region) and the tumor stroma area according to the expression of panCK and DAPI. To obtain multispectral images, the stained slides were scanned using the Mantra System (PerkinElmer, Waltham, Massachusetts, USA), which captures the fluorescent spectra at 20 nm wavelength intervals from 420 to 720 nm with identical exposure time; the scans were combined to build a single stack image. Images of unstained and single-stained sections were used to extract the spectrum of autofluorescence of tissues and each fluorescein, respectively. The extracted images were further used to establish a spectral library required for multispectral unmixing by inForm image analysis software (PerkinElmer, Waltham, Massachusetts, USA). We obtained reconstructed images of sections with the autofluorescence removed by using this spectral library.

Single-cell RNA sequencing and analysis

Paired PL and LM fresh tumor samples were collected from a patient with LM-NSCLC for single-cell RNA sequencing (scRNA-seq). The technological details of single-cell dissociation, scRNA-seq, and analysis were shown in online supplemental materials.

Supplemental material

Statistical analysis

Statistical analyses were performed using the Statistical Package for Social Science (SPSS) software (V.23) and GraphPad Prism software (V.8). The Kaplan-Meier method was used to analyze the survival probability, and the log-rank test was used to calculate the significant differences. The Cox proportional hazard model was applied for the univariate and multivariate analyses to calculate the HRs and 95% CIs. The differential tumor response rate between different organ systems or groups was compared using the χ2 test. Continuous variables were assessed for normality first. A student’s paired t-test was used to analyze the difference between the two groups for normally distributed data, and the Wilcoxon matched-pairs signed-rank test was used for non-normally distributed data. Two-sided p values<0.05 were considered statistically significant.

Results

Response to immunotherapy differs between PL and LM in patients with advanced NSCLC

A total of 52 patients with advanced NSCLC with LM who met the criteria of the clinical cohort were enrolled to evaluate organ-specific tumor response to immunotherapy between PL and LM. Baseline characteristics are summarized in table 1. The median age of patients was 60 (25 to 77), and two-thirds were men. Among them, 67.3% (35/52) were diagnosed with non-squamous NSCLC, and 71.1% of patients (37/52) were treated with ICI as a first-line or second-line treatment strategy. There were 27 patients (51.9%) who received ICIs monotherapy, while the others received ICIs-combination therapy.

Table 1

Patient characteristics (clinical cohort)

All of them had evaluable target lesions in both the liver and lung at baseline, and subsequent imaging assessments were conducted to evaluate organ-specific tumor responses. Figure 1A presents the best percentage change in tumor burden relative to that at baseline in PL and LM, respectively. Among these patients, the OSRR and OSDCR of LM were lower than that of PL (OSRR: 21.1% vs 32.7%, p=0.185; OSDCR: 67.3% vs 86.5%, p=0.020) (figure 1B). Given the influence of ICIs treatment strategies on efficacy in patients with NSCLC, we investigated the differences in the organ-specific tumor response made by different treatment strategies. It seems that the OSRR and OSDCR of LM were numerically lower than that of PL in both the monotherapy group (OSRR: 18.6% vs 29.6%; OSDCR: 59.3% vs 81.5%) and the ICIs-based combination therapy group (OSRR:24.0% vs 36.0%; OSDCR: 76.0% vs 92.0%) (online supplemental figure S1A).

Figure 1

Organ-specific response to immunotherapy differs between PL and LM in patients with advanced NSCLC. (A) Best changes from baseline in target lesions of PL and LM in per patients. (B) Organ-specific response rate and disease control rate in PL and LM. (C) Swimming plot showing progression-free survival and treatment strategy. Pathological type, PD-L1 expression level, gender, and EGFR/ALK mutation status are indicated for each patient. (D) The constituent ratio of different response patterns. (F) Kaplan-Meier estimates of progression-free survival stratified by organ-specific response patterns. ALK, anaplastic lymphoma kinase; ECOG PS, Eastern Cooperative Oncology Group performance status; EGFR, epidermal growth factor receptor; LM, liver metastases; NSCLC, non-small cell lung cancer; OSDCR, organ-specific disease control rate; OSRR, organ-specific response rate; PD, progressive diseases; PD-L1, programmed death-ligand 1; PL, pulmonary lesions; PR, partial response; SD, stable diseases.

Noteworthily, of the total 52 patients, 17 showed mixed response, among which 88.2% (15/17) presented with LM increasing, while the others had the opposite pattern (p=0.002), as shown in figure 1C,D. Univariate and multivariate analyses identified that response pattern is an important feature that was correlated with PFS in addition to the other clinicopathologic variables including performance status, treatment line, and treatment strategies that were correlated with treatment efficacy (online supplemental table S1). In addition, we observed that the presence of mixed response was associated with shorter PFS (median PFS: 2.8 months) compared with those who had disease reduction in both PL and LM (median PFS: 8.9 months, p<0.001) (figure 1E). As for OS, we found that mixed responses presented with inferior OS compared with both reduction groups (median PFS: 9.4 vs 15.1 months, p=0.029). Univariate and multivariate analyses also identified response patterns as an important feature correlated with OS (online supplemental figure S1B and table S2).

Genomic landscapes profiling by WES revealed limited heterogeneity between paired PL and LM

As the responses to immunotherapy differ significantly between PL and LM, we aimed to explore potential cellular and molecular mechanisms underlying the organ-specific response. Therefore, we performed the following translational part by multiomics analysis to obtain genomic, transcriptomic, and immune marker profiles of the 27 pairs of matched PL and LM tumor samples from another independent cohort (online supplemental table S3) following the workflow shown in figure 2A.

Figure 2

Genomic landscapes profiling in paired PL and LM. (A) The workflow chart and samples grouping of the translational part. (B) Venn diagram of genetic variation, SNV, and CNV in the PL and LM groups. (C) The comparison of CNV all, CNV loss, and CNV gain burden between the PL and LM groups. (D) The comparison of TMB between the PL and LM groups. (E) The comparison of neoantigens number predicted by SNV and indel(-type) variants between the PL and LM groups. (F) The most common mutational genes are depicted in the heat map. Note: *, p<0.05. CNV, copy number variation; IHC, immunohistochemistry; LM, liver metastases; LUAD, lung adenocarcinoma; LUSC, lung squamous carcinoma; mIHC, multiplex IHC; muts/Mb, mutations per megabase; NSCLC, non-small cell lung cancer; PL, pulmonary lesions; SNV, single nucleotide variant; TMB, tumor mutational burden; WES, whole-exome sequencing.

To explore the mutational heterogeneity, whole-exome sequencing (WES) was performed to identify somatic mutations in six pairs of tumor samples (including three PL−LM and three PL=LM to ensure the paired samples are in the same systemic treatment status as far as possible). We first compared the consistency of somatic mutations between PL and LM. The consistency of genetic variation detecting by WES was 47.5%, among which SNV was 41.1% and copy number variation (CNV) was 52.8% (figure 2B). Although the CNV loss burden was slightly higher in LM than in PL, there was no significant difference between these two groups (figure 2C).

Regarding the tumor mutational burden (TMB), the mean value of TMB of all the lesions was 6.1 mutations/Mb, ranging from 1.4 to 11.0 mutations/Mb in PL and from 3.4 to 12.5 mutations/Mb in LM (figure 2D). No significant difference was found in TMB, and the number of neoantigens predicted by SNV and indel(-type) variants between PL and LM (figure 2E).

We further analyzed the heterogeneity of mutational genes, finding no significantly different gene mutations between PL and LM. The most common gene mutations occurred in FLG2 (100% in PL and 83.3% in LM), MUC17 (100% in PL and 83.3% in LM), and KLF18 (100% in PL and 66.7% in LM), which showed a high degree of agreement between the two groups. Comparison of PL and LM in variant classification and SNV class revealed high consistency in both groups. Either in PL or LM, the most common variant classification was missense mutation, and the SNV class was C>T (figure 2F, online supplemental figure S2).

Pro-immune signaling pathways were enriched in PL than LM by bulk-transcriptome sequencing

To evaluate transcriptional heterogeneity between PL and LM, we performed RNA-sequencing in 12 matched tumor samples collected simultaneously from six patients with NSCLC. We first compared the gene expression profiles, identifying 571 differentially expressed genes (DEGs), comprizing 469 genes upregulated in LM and 102 in PL (figure 3A), which was indicative of transcriptional intertumoral heterogeneity in NSCLC within individual tumors. By further analyzing all downregulated and upregulated DEGs from each group, we found that genes overexpressed in LM were mainly related to metabolic enzymes, which was consistent with the normal physiological function of the liver. In contrast, some immune-related genes, including CD8A, LCK, and ICOS, which were relevant to the activation of T cells, were downregulated in LM.

Figure 3

LM shows an immunosuppressive tumor microenvironment. (A) A volcano plot of differentially expressed genes between PL and LM with four immune-related genes downregulated in LM marked in a frame. (B) Gene set enrichment analysis between PL and LM among three immunotherapy efficacy-related pathways: the antigen processing and presentation, NK cell-mediated cytotoxicity, and T-cell receptor signaling pathway. (C) Gene Ontology analysis of differentially enriched biological processes upregulated (red) or downregulated (blue) in LM. (D) The comparison of 22 immune cells proportion estimated by the CIBERSORT approach between the PL and LM groups. Note: *, p<0.05. DC, dendritic cell; KEGG, Kyoto encyclopedia of genes and genomes; LM, liver metastases; NK, natural killer; PL, pulmonary lesions.

Several gene signatures that could influence the efficacy of immunotherapy were examined. Three gene signatures positively associated with immunotherapy efficacy including the antigen processing and presentation, NK cell-mediated cytotoxicity, and T-cell receptor signaling pathway were slightly enriched in PL, but not LM (figure 3B). Gene ontology analysis of differentially enriched biological processes also showed LM was significantly involved in the metabolic and catabolic processes of various substances. In contrast, PL was significantly involved in immune cell activation and immune response process compared with LM (figure 3C), which could partially explain the different responses to immunotherapy.

Next, we used the CIBERSORT approach to preliminarily explore the patterns of immune cell infiltration. We found that of the 22 immune cell subtypes involved in the CIBERSORT approach, only CD8+T cells infiltration significantly differed between these two groups (figure 3D). Using gene expression data to estimate the abundances of member cell types in a mixed cell population, we found that more CD8+T cells were infiltered in PL rather than LM (p=0.022). Moreover, an interesting trend was observed where the infiltration of activated memory CD4+T cell was higher in the PL group than in LM (p=0.082), while the infiltration of the resting memory CD4+T cell was reversed (p=0.059).

Lower fractions of T cells and NK cells were infiltrated in LM than in PL by mIHC

To further confirm the patterns of immune cell infiltration and explore the spatial distribution of the immune microenvironment between the PL and LM at the protein level, the mIHC was conducted in 24 pairs of tumor samples. As it was reported that systemic anticancer therapies have a dramatic impact on TME, we included four PT−LM and four PT=LM in a subgroup analysis, named group A, to exclude the influence of systemic treatment as much as possible. Herein, five morphology markers including panCK, CD8A, CD56, CD68, and HLA-DR were used to identify panCK+ tumor cells, CD8+ T cells, CD56dim+ cytotoxic NK cells, CD56bright+ immunoregulatory NK cells, CD68+HLA-DR+ macrophages (M1-like), and CD68+HLA-DR macrophages (M2-like).

Comparison of the infiltration of immune cells between the parenchyma and stroma area revealed higher intratumor heterogeneity of immune cell distribution in LM than in PL, and such heterogeneity was more obvious in group A. By analyzing the relative proportion separately, we found the infiltration trend of CD8+, CD56+, and CD68+ cells were similar in the parenchyma area and stroma area in the PL group, while CD56+ cells occupied most of the LM-tumor parenchyma area and CD68+ cells of LM-tumor stroma area (figure 4A,B).

Figure 4

Lower fractions of CD8+ cell and CD56dim+ cell infiltrated in LM. (A) The positive rate of CD8+ cells, CD56dim+ cells, CD56bright+ cells, CD68+HLA-DR+ cells, and CD68+HLA-DR cells in the total tumor area, tumor parenchyma area and tumor stroma area in the PL and LM groups. (B) The relative infiltration proportion of CD8+ cells, CD56dim+ cells, CD56bright+ cells, CD68+HLA-DR+ cells, and CD68+HLA-DR cells in the tumor parenchyma area and tumor stroma area in the PL and LM groups. (C) The changes of CD8+ cells, CD56dim+ cells, CD56bright+ cells, CD68+HLA-DR+ cells, and CD68+HLA-DR cells infiltration from PL to LM in the total population (24 pairs) and group A (8 pairs). (D) The changes of CD8+ cells and CD56dim+ cells in each patient of group A. Note: *, p<0.05. HLA-DR, human leukocyte antigen-DR; LM, liver metastases; PL, pulmonary lesions.

Next, we analyzed the differences between PL and LM. It was found that, in the total population (n=24 pairs), PL and LM share similar patterns of the relative distribution of immune cells in whether parenchymal or stroma area, for example, CD68+HLA-DR macrophages occupy most of the stroma area in both the PL and LM group (table 2). In group A (n=8 pairs), a more homogeneous group, higher heterogeneity was observed between these groups (figure 4A). It illustrated that the effect of systemic treatment on the TME is substantial, which may mask the actual differences between the PL and LM groups. Similarly, when we compared the density and positive rate of these immune markers between these two groups in the total population, only CD56+ cells were significantly higher in the PL group compared with the LM group. No significant difference was found in CD8+, CD68+HLA-DR+, and CD68+HLA-DR cells (table 2). Subgroup analysis of group A showed significantly lower fractions of CD8+ cell and CD56dim+cell in LM, which was consistent with our findings extrapolated from the bulk-transcriptome sequencing data (figure 4C,D, online supplemental figure S3,S4).

Table 2

Comparison of the expression of multiple immune markers between PL and LM (total population)

To examine the correlation between immune cell infiltration and response to immunotherapy, we analyzed the treatment history of patients in the translational part. Seven patients received immunotherapy with recorded imaging efficacy evaluation, among whom three showed mixed responses to immunotherapy. All these three patients presented with PL reduction while LM increasing, and exhibited lower infiltration of CD8+ cells in LM compared with PL (online supplemental table S4).

Higher infiltration of effector CD8+ T cells and cytotoxicity of NK cells in PL by scRNA-seq

Then we performed scRNA-seq in paired PL and LM tumor samples from an adult patient with NSCLC. Fresh tumor samples were both collected at progression to third epidermal growth factor receptor-tyrosine kinase inhibitor (EGFR)-TKI (online supplemental figure S5A). As we found that the LM TME displays lower levels of T and NK cell activation and infiltration than PL from bulk-RNA sequencing and mIHC, the comparison of gene expression in T/NK cells at the single-cell level focused on further validating these previous findings.

A total of 22,374 cells were identified, among which 5092 were annotated as T/NK cells according to the expression of canonical markers (CD3D, CD2, TRAC, TRBC2). Reclustering of NK/T cells revealed eight clusters, which were designated as C1-effector CD8 T cells (CD8Teff), C2-mucosal associated invariant CD8 T cells (CD8MAIT), C3-effector memory CD8 T cells (CD8Tem), C4-NK cells, C5-Helper CD4 T cells (HelperT), C6-Regulatory CD4 T cells (Treg), C7-ISG T cells (IFN-associated-T, highly expressed interferon-stimulated genes), and C8-MKI67 T cells (proliferatingT, highly expressed cell-cycle genes) (figure 5A,B). Consistent with the previous findings, the proportions of Treg and CD8Tem were higher in LM, while CD8Teff and HelperT were higher in PL (figure 5C,D).

Figure 5

The scRNA-seq profiling of T/NK cells in paired PL and LM. (A) Uniform manifold approximation and projection (UMAP) plot of T/NK cells profiled in the present study colored by subcluster. (B) Feature plots showing the normalized expression of canonical marker genes in each subcluster of T/NK cells. (C) Proportions of T/NK cell subtypes in PL and LM, respectively. (D) UMAP plot of the representative gene for CD8Teff, HelperT, Treg, and CD8Tem. (E) GSEA between PL and LM among the T-cell receptor signaling pathway and NK cell-mediated cytotoxicity in isolated T and NK cells. (F) UMAP plot showing the T-cell activation and cytotoxicity scores in PL and LM. (G) Column chart showing the T-cell activation and cytotoxicity scores in each subcluster of T/NK cells. (H) Bubble diagram depicting the significant chemokines and (I) ligand–receptor interactions in PL and LM. GSEA, gene set enrichment analysis; LM, liver metastases; NK, natural killer; PL, pulmonary lesions; scRNA-seq, single-cell RNA sequencing.

NK cell-mediated cytotoxicity and T-cell receptor signaling pathway were also enriched in PL compared with LM when NK cells and T cells were isolated for analysis at the single-cell level (figure 5E), after which T-cell activation and cytotoxicity scores were evaluated for each T/NK cluster. All eight clusters exhibited higher activation scores in PL than in LM, among which CD8Teff and IFN-associated-T showed the highest score. NK cells, HelperT, IFN-associated-T, and proliferatingT possessed higher cytotoxicity scores in PL than in LM. Noteworthily, NK cells had substantially higher cytotoxicity scores than the other clusters (figure 5F,G).

Understanding the molecular interactions is important to comprehend the mechanisms of distinct T/NK cells distribution. We then performed the cell–cell connections in PL and LM, respectively, finding that PL has higher levels of Teff-recruiting chemokines (CCL23, CCL13, CCL18, and CCL3), while LM could recruit more Treg through CCL18_CCR8 and SPP1_CCR8 (figure 5H).18 Besides, T/NK cells in LM could be regulated through specific ligand–receptor pairs: such as CLEC2D_ KLRB1 and SIGLEC10_CD52, which were both associated with the inhibition of T/NK cells proliferation and activation.19 On the other hand, the PDCD1 and CD274 pair, which was the target of PD-1/L1 inhibitors, was detected among the interactions of CD8Teff with macrophages in PL instead of LM (figure 5I). It was suspected that the activity of CD8Teff in LM-TME could not be restored after the usage of PD-1/L1 inhibitors. Besides, we found that programmed cell death pathways in T cells, including apoptosis, ferroptosis, and necroptosis, were all enriched in LM instead of PL (online supplemental figure S5).

Overall, our results collectively illustrated a more suppressive TME in LM compared with PL: lower levels of activation and infiltration of CD8+T and NK cells, and diminution of their cytotoxic capacity.

Discussion

It is the first study using strictly paired tumor samples to explore organ-specific response and TME between PL and LM from patients with NSCLC, which revealed organ-specific differences in the clinical tumor response to ICIs and its protentional mechanism across genomic, transcriptomic, and multiple immune markers protein levels. We demonstrated limited heterogeneity in the genomic landscape between paired PL and LM tumors; however, the LM TME displayed lower immune activation and infiltration levels than PL. The more suppressive TME profiling by multiomics analysis might explain the inferior organ-specific response of immunotherapy in LM.

Previous clinical studies on NSCLC showed that the presence of baseline LM was associated with diminished efficacy of immunotherapy.20 21 Our study included a clinical cohort to evaluate organ-specific tumor responses, finding that the responses to ICIs differed significantly between PL and LM. An inferior response to ICIs in LM and a favorable tumor response in PL was observed. Consistent with our findings, some studies reported that tumor responses to ICIs varied significantly depending on the organ in which the tumor was located, which is not unique to NSCLC.10 21 In a subgroup analysis from KEYNOTE-001, individual lesions from 37 patients with melanoma were analyzed, showing that lung lesions had the highest rate of CR (42.3%), followed by peritoneal (37.3%) and liver (24.4%) metastatic lesions.9 Moreover, we found that such a mixed organ-specific response to ICIs was significantly correlated with negative survival benefits.

Mixed response is not rare in patients with NSCLC treated with TKIs or chemotherapy. It was reported that the incidence of a mixed response to EGFR-TKI and chemotherapy was 21.5% (53/246), which mainly originated from genetic heterogeneity.22 The present study found genomic landscapes profiling by WES revealing limited heterogeneity between paired PL and LM tumors. High consistency in SNV, CNV, TMB, mutational genes, and variant classification was also detected. Accordingly, we concluded that genetic heterogeneity could not explain the different organ-specific responses of immunotherapy between paired PL and LM.

TME plays an important role in ICIs inducing antitumor effects. Heterogeneous TME of various organs may potentially influence tumor growth and responses to anticancer treatment.6 Our findings further clarify the tissue specificity of TME in NSCLC with LM. Various cancer models found differences in the tumor and TME were also metastatic organ dependent.11 23–25 It was reported that tumors from renal, colon, or prostate cell lines in orthotopic locations responded to immunotherapy much less than the same tumor type subcutaneously located despite the injection of genetically-matched cancer cells. Besides, orthotopic tumors had a microenvironment related to immunosuppression compared with subcutaneous tumors. It was emphasized that the site of disease determines the tumor response to immunotherapy by cross-implanting renal tumor cells harvested from one site into another.26 All these findings indicate that the tissue surrounding the tumor, varying with the anatomical site, has a major impact on the TME.

Further investigations are needed to identify distinctive tissue-specific signatures important in regulating the TME, as these could potentially inform future therapeutic strategies against hepatic-metastatic NSCLC. However, studying these organ-specific differences has been challenging, for collecting matched tumor samples from different metastatic sites of the same individual is inherently tricky. Since there is rarely ever a need for simultaneous multisite biopsy for clinical decision-making in NSCLC. Previous studies about organ-specific TME have always used animal models or non-strictly paired NSCLC tumor samples.11 12 Comparative analysis of paired tumor samples from the same patients could be more indicative of the organ-specific TME differences after excluding individual heterogeneity. Herein, we detected 27 pairs of strictly paired PL and LM tumor samples, including six simultaneously collected pairs, from the same individuals to identify the organ-specific TME signatures of PL and LM in NSCLC.

Because of its inhibitory immunomodulation properties, the liver is viewed as a unique organ constantly exposed to various foreign antigens released from the gastrointestinal tract.27 28 Previous studies on LM TME were mostly in colorectal cancer since patients with colorectal cancer with LM still have the opportunity for radical surgical resection against the primary lesion and LM. However, data about NSCLC-LM is lacking, especially those comparing LM and PL. In the present study, we suggested a higher threshold of immune suppression to overcome in LM for antitumor immunotherapy, including lower levels of activation and infiltration of CD8+T cells as well as NK cells and diminution of their cytotoxic capacity in the overall comparison of TME in paired PL and LM. When examining the correlation between immune cell infiltration and response to immunotherapy, we also found that those with mixed response infiltrated a lower fraction of CD8+ cells in LM than in PL. It was suspected that CD8+ T cells might have an important role in the inferior organ-specific response of LM. Adoptive T or NK cell therapy might be a potentially effective approach to improve the infiltration and activity of T cells and NK cells in the LM. A combination of ICIs and adoptive T/NK cell therapies are currently being investigated in clinical trials (NCT03645928, NCT03215810, NCT05395052, NCT03941262), and some of them had preliminarily reported encouraging results.29 30

In addition, CD274 and PDCD1_CD274 interaction expressions were significantly lower in LM compared with PL. It is possible that PD-1/L1 blockade as a single agent might not exert a satisfactory local antitumor effect in LM, and a combination of anti-PD1 therapy with strategies to improve the TME might be a better solution. Growing evidence has revealed that angiogenic molecules, particularly vascular endothelial growth factor, have an important role in immunosuppression in the TME via various mechanisms.31 Several studies reported that anti-angiogenic therapy could increase the infiltration of cytotoxic T lymphocytes and enhance the expression of PD-L1.32 33 A randomized trial comparing the efficacy of PD-L1 inhibitors with or without anti-angiogenesis drugs targeting patients with NSCLC with LM is currently recruiting (NCT04563338).

However, the present study has some limitations. Given the nature of a single-center observational clinical cohort study, our results must be interpreted cautiously. In addition, the PL and LM samples we chose in this study were strictly paired samples from the same individuals, which resulted in a limited sample size (n=27 pairs, including 6 simultaneously collected pairs). Still, comparative analysis of strictly paired tumor samples could be more indicative of the organ-specific TME differences after excluding individual heterogeneity. Even limited, it is the first study using a relatively large sample size to explore organ-specific TME between PL and LM by analyzing paired tumor samples from patients with NSCLC. Multiomic approaches, including WES, RNA-sequencing, and mIHC were used to extract as much information about TME as possible from limited samples; however, more studies are needed to confirm our results further.

In summary, our study demonstrated the organ-specific tumor response to ICIs and delineated tissue specificity of TME in pulmonary and liver-metastatic NSCLC. A more immunosuppressive tumor ecosystem of LM in NSCLC was uncovered, which provides insights into the mechanisms underlying the inferior organ-specific response of immunotherapy and information for the development of precise immunotherapy strategies for patients with NSCLC with LM.

Data availability statement

Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information.

Ethics statements

Patient consent for publication

Ethics approval

This study involves human participants and was approved by the Research Ethics Committee of Guangdong Provincial People’s Hospital (No. GDREC2016262H). Participants gave informed consent to participate in the study before taking part.

Acknowledgments

The authors would like to thank the patients, their families and the study personnel who participated in this study. The authors thank the Nanjing Tongyuan Medical Laboratory for their help with the single-cell RNA sequencing analysis.

References

Supplementary materials

  • Supplementary Data

    This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.

Footnotes

  • Twitter @xucr

  • J-YD and QG contributed equally.

  • Contributors Conception and design: QZ and J-YD. Acquisition tumor samples of patients in Guangdong Lung Cancer Institute: J-YD, QG, M-YY, X-RY, XW, JS, Z-HC and WZ. Development of methodology: J-YD, QG, LY and C-RX. Analysis and interpretation of data: J-YD, QG, LY and H-HY. Writing and revision of the manuscript: J-YD, QG and QZ. Review of the manuscript and put forward suggestions for amendment: All authors. Final approval of the version to be submitted: All authors. QZ is responsible for the overall content as guarantor.

  • Funding The study was supported by the National Natural Science Foundation of China (Grant No. 82072562 to QZ) and the High-level Hospital Construction Project (Grant No. DFJH201810 to QZ).

  • Competing interests QZ reports honoraria from AstraZeneca, Boehringer Ingelheim, BMS, EliLilly, MSD, Pfizer, Roche, and Sanofi, outside the submitted work. Y-LW reports personal financial interests: consulting and advisory services, speaking engagements of Roche, AstraZeneca, EliLilly, Boehringer Ingelheim, Sanofi, MSD, and BMS. LY is an employee of Nanjing Geenseeq Technology. All remaining authors declared no competing interests.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.