Article Text

Original research
High-throughput transcriptome profiling indicates ribosomal RNAs to be associated with resistance to immunotherapy in non-small cell lung cancer (NSCLC)
  1. Myrto K Moutafi1,
  2. Katherine M Bates1,
  3. Thazin Nwe Aung2,3,
  4. Rolando Garcia Milian4,
  5. Vasiliki Xirou2,3,
  6. Ioannis A Vathiotis2,
  7. Niki Gavrielatou2,3,
  8. Athanasios Angelakis5,6,
  9. Kurt A Schalper2,
  10. Leonidas Salichos7 and
  11. David L Rimm1
  1. 1Department of Pathology, Yale University School of Medicine, New Haven, Connecticut, USA
  2. 2Department of Pathology, Yale School of Medicine, New Haven, Connecticut, USA
  3. 3Yale School of Medicine, New Haven, Connecticut, USA
  4. 4Bioinformatics Support Program, Cushing/Whitney Medical Library, Yale School of Medicine, New Haven, Connecticut, USA
  5. 5Epidemiology and Data Science, Amsterdam UMC Locatie AMC, Amsterdam, The Netherlands
  6. 6Department of Methodology, Amsterdam Public Health Research Institute, Amsterdam, The Netherlands
  7. 7Biomedical Data Science Center Director, Center for Cancer Research, Department of Computational Biology at New York Institute of Technology, New York Institute of Technology, Old Westbury, New York, USA
  1. Correspondence to Dr David L Rimm; david.rimm{at}yale.edu

Abstract

Background Despite the impressive outcomes with immune checkpoint inhibitor (ICI) in non-small cell lung cancer (NSCLC), only a minority of the patients show long-term benefits from ICI. In this study, we used retrospective cohorts of ICI treated patients with NSCLC to discover and validate spatially resolved protein markers associated with resistance to programmed cell death protein-1 (PD-1) axis inhibition.

Methods Pretreatment samples from 56 patients with NSCLC treated with ICI were collected and analyzed in a tissue microarray (TMA) format in including four different tumor regions per patient using the GeoMx platform for spatially informed transcriptomics. 34 patients had assessable tissue with tumor compartment in all 4 TMA spots, 22 with leukocyte compartment and 12 with CD68 compartment. The patients’ tissue that was not assessable in fourfold redundancy in each compartment was designated as the validation cohort; cytokeratin (CK) (N=22), leukocytes CD45 (N=31), macrophages, CD68 (N=43). The human whole transcriptome, represented by~18,000 individual genes assessed by oligonucleotide-tagged in situ hybridization, was sequenced on the NovaSeq platform to quantify the RNAs present in each region of interest.

Results 54,000 gene variables were generated per case, from them 25,740 were analyzed after removing targets with expression lower than a prespecified frequency. Cox proportional-hazards model analysis was performed for overall and progression-free survival (OS, PFS, respectively). After identifying genes significantly associated with limited survival benefit (HR>1)/progression per spot per patient, we used the intersection of them across the four TMA spots per patient. This resulted in a list of 12 genes in the tumor-cell compartment (RPL13A, GNL3, FAM83A, CYBA, ACSL4, SLC25A6, EPAS1, RPL5, APOL1, HSPD1, RPS4Y1, ADI1). RPL13A, GNL3 in tumor-cell compartment were also significantly associated with OS and PFS, respectively, in the validation cohort (CK: HR, 2.48; p=0.02 and HR, 5.33; p=0.04). In CD45 compartment, secreted frizzled-related protein 2, was associated with OS in the discovery cohort but not in the validation cohort. Similarly, in the CD68 compartment ARHGAP and PNN interacting serine and arginine rich protein were significantly associated with PFS and OS, respectively, in the majority but not all four spots per patient.

Conclusion This work highlights RPL13A and GNL3 as potential indicative biomarkers of resistance to PD-1 axis blockade that might help to improve precision immunotherapy strategies for lung cancer.

  • immune checkpoint inhibitor
  • lung cancer
  • biomarker

Data availability statement

Data are available upon reasonable request.

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

  • Inherent or acquired resistance to programmed cell death protein-1/programmed death-ligand 1 blockade treatment results from the persistence or emergence of tumor cell subsets often with metastatic capabilities.

WHAT THIS STUDY ADDS

  • Using molecular compartment-specific spatial transcriptomics on four cores from different tumor areas of the same patient tumor sample, we have discovered two candidate, tumor-compartment (cytokeratin+ positive cells) specific genes for immune checkpoint inhibitor resistance (RPL13A and GNL3).

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • Results from this study may lead to a novel, spatially defined, transcriptomic approach for developing new biomarkers for immunotherapy resistance.

Introduction

Immunotherapies that block the programmed cell death protein-1/programmed death-ligand 1 (PD-1/PD-L1) checkpoint have transformed lung cancer treatment and outcomes.1 However, most patients do not derive significant clinical benefits and patient selection biomarkers for immune checkpoint inhibitors (ICI) are needed.2 The mechanisms underlying such variable responses remain incompletely understood. A better understanding of such mechanisms could be used to broaden the spectrum of responding patients and enhance the rate of ICI response. Therefore, there is a medical need to identify biomarkers or mechanisms by which tumors show resistance to PD-1 axis inhibition, to expand the treatment opportunities and to select patients who may opt out of immunotherapy treatment in the adjuvant setting.

Recently, ICI have been approved in the perioperative setting in melanoma and non-small cell lung cancer (NSCLC).3 4 This presents a new challenge for oncologists. Since the ICIs are associated with a series of side effects, the oncologist must decide if prescribing an ICI to a patient may already be cured by surgery alone. For example, in NSCLC there is a 10% increase in disease-free survival for those given atezolizumab in the IMpower010 trial, compared with the placebo arm.5 6 However, there was a high rate of thyroiditis and a 1.5% rate of death in this population. Similar results have been seen using ICIs alone or in combination with chemotherapy in the neoadjuvant setting.7 Thus, the oncologist and patient must make a difficult decision regarding the benefit of perioperative therapy. This decision would be substantially assisted by biomarker testing that would estimate either the prognosis of the disease or, possibly more easily, the likelihood of resistance to immunotherapy. If the biomarker showed that the specific tumor was unlikely to benefit from ICI, then the patient might decide to forego treatment and the attendant risk of adverse events.

Finding such a biomarker of ICI resistance has been a challenge. PD-L1 immunohistochemistry (IHC) assays are neither sensitive nor specific with an area under the curve (AUC) less than 0.7.8 9 Multigene immune signatures and tumor mutation burden have been tested to address the relatively low PD-L1 IHC AUC.10–14 However, these biomarkers show variable accuracy and their applicability in clinical practice is not established.15 In the PIONeeR trial, a comprehensive multiparametric and longitudinal biomarker assessment of tumor tissue and blood was used to assess response to first-line pembrolizumab plus platinum-based chemotherapy. The signature showed a high predictive performance for progression-free survival (PFS) at 1 year on a training set of data (n=116) (C-index 0.77 and AUC 0.80) but achieved lower performance in a testing data set (n=39) (C-index 0.73 and AUC 0.68) (ESMO IO2022, Abstract 3MO).16 Predictive models of response/resistance to ICI in advanced NSCLC still need further validation to tailor treatment strategies.

Part of the reason for the absence of effective biomarkers of resistance may be the complex nature of tumors and the inability of conventional methods to define molecular compartment-specific biomarkers. The different role of tumor and immune cell populations can be revealed by performing a spatially informed analysis of biomarkers for ICI.17 Here we performed a complete spatial transcriptomic study attempting to discover new RNA biomarkers for resistance to ICI therapy using a spatially informed high-plex discovery tool, the GeoMx Digital Spatial Profiling (DSP) System. This tool provides an unsupervised measurement of 18,000+protein-coding genes and enables their assessment in a spatially informed manner in single formalin-fixed paraffin embedded (FFPE) tissue sections by using next-generation sequencing (NGS). The DSP platform measures the quantity of RNA by counting unique indexing oligonucleotides assigned to each target of interest. Areas of interest (AOIs) or molecularly defined spatial compartments of any shape, including non-contiguous areas, can be profiled in a region-specific manner.18 19 In this study, we describe the use of DSP technology as a discovery tool to find spatially resolved RNA markers associated with resistance to ICI in advanced NSCLC in three molecularly defined spatial compartments (tumor-cytokeratin, CK+, leukocytes-CD45+/CD68−, and macrophages-CD45+/CD68+). Then, among candidate predictors per compartment, we validated the association with outcome using independent cohorts.

Materials and methods

Patient cohorts and tissue microarrays

We analyzed retrospectively collected baseline/pretreatment FFPE tumor specimens represented in four tissue microarrays (TMA) format from 56 patients with NSCLC treated with ICI in the advanced setting between 2017 and 2019 at Yale University School of Medicine, New Haven, Connecticut, USA (Yale Tissue MicroArray 471; YTMA471). Tumor samples were collected from retrospective immunotherapy-treated cohorts from a single institution, not from randomized clinical trials. Tissue samples from these 56 patients were used to create two independent sets; a discovery cohort and a validation set, for each molecular compartment, for the discovery of biomarkers of resistance to immunotherapy (table 1). All tissue samples were collected and used under the approval of the Yale Human Research Protection Program (HRPP) Institutional Review Board). The Yale HRPP approved the patient informed consent or in some cases waiver of consent all in accordance with the ethical guidelines of the US Common Rule.

Table 1

Clinicopathologic characteristics of the cohorts analyzed

For TMA construction, tumors were reviewed by a pathologist using H&E–stained preparations to select representative tumor areas. Then, four cores (area; 0.28 mm2, diameter; 600 µm each) were extracted from the region selected by the pathologist from each tumor block and arrayed in four recipient TMA master blocks (TMA 1–4). Each TMA block thus containing one 0.6 mm diameter tumor core per NSCLC case (figure 1, online supplemental figure 1). Tumor core selection was based on the presence of viable tumor (CK+) cells with morphology and associated stromal features representing the disease as judged by a board-certified pathology. For all the experiments, we assessed four slides derived from four independently constructed YTMA471 blocks, each block containing one tumor core per patient.

Supplemental material

Supplemental material

Supplemental material

Figure 1

Workflow. (A) Four independent YTMA471 blocks, were analyzed; each block containing one non-adjacent tumor core per patient. (B) For the CK compartment 34 cores had evaluable CK+tissue in all four YTMA blocks (marked with magenta); 22 cores had evaluable CK+ tissue in at least one block and consisted the validation set. (C) Four independent Cox survival analysis were performed for each block. The intersection of the genes significantly associated with survival is shown in the middle. (D) The intersection of the 12 genes of the fourfold redundant cohorts with the validation set for the CK compartment consisted of 2 genes. (E) Kaplan-Meier graph of one of the two genes in the validation cohort. CK, cytokeratin; NSCLC, non-small cell lung cancer.

Figure 2

Volcano plots. (A) Volcano plot showing genes associations with survival in Block 1, (B) Block 2, (C) Block 3, (D) Block 4× axis logHR; HR>1 indicates shorter survival. (E) Taking the intersection of all the significantly associated RNAs with overall survival (OS) throughout the four YTMAs, RPL13A, GNL3, FAM83A, CYBA, ACSL4, SLC25A6, EPAS1, RPL5, APOL1, HSPD1, RPS4Y1 and ADI1 were correlated with shorter OS regardless the tumor core (HR>1, PLogRank <0.05). (F) Two of them ribosomal protein L13a (RPL13A), GNL3 (both participating in ribosomal RNA processing in the nucleolus and cytosol) they were also associated with shorter OS in the validation cohort (HR 2.48, PLogRank 0.02; HR 5.33, PLogRank 0.04). CK, cytokeratin.

Figure 3

Shown are Kaplan-Meier estimates of overall and progression-free survival, for CK compartment in the validation cohort. (A) Overall survival—RPL13A, (B). Overall survival—GNL3, (C). Progression-fee survival—RPL13A, (D) progression-free survival—GNL3. CK, cytokeratin; OS, overall survival; PFS, progression-fee survival.

Our initial cohort included 56 NSCLC cases, each represented by four TMA cores in four master blocks. For each compartment, two sets were created. The discovery set included patient samples that were available in all four cores. The validation set included the rest of the patient samples. For example, for CK compartment samples from 34 patients (out of 56) were evaluable and these cases defined the CK-discovery cohort. One evaluable core per case derived from the rest of the 22 patient samples defined the validation set. For CD45 compartment samples from 22 patients (out of 56) were evaluable in all four cores and these cases defined the CD45-discovery cohort. One evaluable core per case derived from the rest 31 patient samples defined the validation set. For the CD68 compartment, samples from 12 patients (out of 56) were evaluable in all four cores and these cases defined the CD68-discovery cohort. One evaluable core per block derived from the rest 43 patient samples defined the validation set (table 1). Demographic characteristics were similarly distributed between the two cohorts. We used Response Evaluation Criteria In Solid Tumors V.1.1 to assess treatment response to immune checkpoint blockade. We defined durable clinical benefit (DCB) as having experienced partial response or stable disease lasting ≥6 months as the best response, whereas non-durable clinical benefit was defined as a progressive disease or stable disease lasting <6 months. Overall survival (OS) and PFS were calculated from the treatment start date to the date of death or loss of follow-up, or the date of disease progression, death or loss of follow-up, respectively.

Digital Spatial Profiling

Briefly, once the slides were deparaffined and subjected to antigen retrieval procedures, we incubated them overnight with three fluorescent-labeled compartment visualization antibodies to detect tumor cells (pan-CK), leukocytes (CD45), and macrophages (CD68) to define the spatially unique, molecular compartments, as previously described.20

These were overlayed with a cocktail of 18,815 unique photocleavable oligonucleotide-conjugated in situ hybridization (ISH) probes targeting more than 18,815 protein-coding genes previously validated by NanoString as previously shown.21 Once the staining was completed, slides were loaded onto the GeoMx DSP instrument, where they were scanned to produce a digital fluorescent image of the tissue.

Individual regions of interest (ROI) of a maximum diameter of 600 µm covering the entire TMA core were generated, and then each ROI was segmented into three molecularly defined tissue compartments by fluorescent colocalization: tumor compartment (CK+), leukocyte compartment (CD45+/ CD68−), and macrophage compartment (CD68+) (fluorescently labeled antibodies to PanCK, CD45, CD68 identify epithelium, immune cells and macrophages, respectively, and a nucleic acid stain identifies nuclei).

Every ROI was sequentially exposed to ultraviolet light of a programmable digital micromirror device, that precisely illuminates the ROIs to decouple the oligonucleotides from the profiling reagent. Photocleaved, decoupled, oligos, each containing a readout tag sequence identifier (RTS ID) specific to a biological target, were collected via microcapillary aspiration and deposited into a 96-well plate. These RTS ID were then analyzed following a pooling, purification, and quality control process, and the final NGS library was sequenced on a NovaSeq platform at Yale Center for Genome Analysis. Sequencing libraries were generated by PCR from the photo-released indexing oligos and AOI-specific Illumina adapter sequences, and unique i5 and i7 sample indices were added. Pooled libraries were paired-end sequenced at 2×27 base pairs and with a unique dual indexes workflow on a NovaSeq S2 Flow cell instrument. Next, the FASTQ files generated from the sequencing run were processed into digital count conversion (DCC) files using the NanoString GeoMx NGS Pipeline. The DCC files were used finally for data analysis in our study.

Messenger RNA in situ hybridization

ISH for GNL3 and RPL13A messenger RNA was performed using the RNAscope Multiplex Fluorescent V.2 kit according to the manufacturer’s instructions, with modifications to allow immunofluorescent staining of CK protein. Briefly, 5 µm thick TMA sections were subjected to heat-mediated antigen retrieval and protease digestion, followed by hybridization with a target probes to GNL3, RPL13A, and the bacterial gene DapB (negative control). Following amplification and HRP binding, the GNL3 and RPL13A hybridization signals were labeled with Cy5-tyramide. After HRP blocking, sections were incubated with a wide-spectrum monoclonal mouse CK antibody (AE1/AE3) directly conjugated to Alexa Fluor 488 for 30 min at 40°C (1:100, Invitrogen, Waltham, Massachusett, USA). Slides were then counterstained with using 4,6-diamidino-2-phenylindole and mounted with ProLong Gold.

The RNAscope assay for GNL3 and RPL13A was performed on serial sections of YTMA471, which was the cohort used for biomarkers discovery, and GNL3 and RPL13A expression was quantified using QYMIA.22 QYMIA is a molecular compartment-based quantification method that measures the target by dividing the sum of the target’s signal intensity within the compartment by the sum of the pixel area of the compartment comparable to the AQUA method of protein quantification on slides. In QuPath, we used wide-spectrum CK fluorescent staining to define the tumor compartment, and then used QYMIA to score GNL3 and RPL13A expression within the tumor area. The resulting QYMIA scores were then correlated with the DSP counts (online supplemental figure S2).

Data analysis pipeline

We first used the GeoMx DSP Control Center to analyze the four data sets created from the four blocks used in this study. The sequence runs were assessed using the same thresholds and settings for all four blocks. Initial data processing and sample-based quality control (QC) were conducted using GeoMx DSP Control Center, any ROIs that were assigned with QC flags indicating low qualities (including “Low Percent Aligned Reads”= “80%”, “Low Percent Stitched Reads”=“80%”, “Low Surface Area”=“16 000 µm2”, “Low Nuclei Count”=“3”) were excluded from further analysis. Negative probes were also removed from the data set in the analyses since they were used in the generation of probe QC count with an outlier-detection strategy.

Data QC was performed, and filter parameters allowed to establish an expression threshold and a frequency for targets to filter out those with an expression lower than a prespecified limit of quantification (≤2 SD above the geometric mean of negative probes) in at least 5% of segments.23 This approach aims to identify genes that are lowly expressed in over 95% of the ROIs; such genes were then removed from the analysis because genes with consistently low expression are unlikely to be determined as significantly differential expressed genes given their low level of expression. Data, then, were scaled using third quantile (Q3) normalization, which is a method using the 75th quantile as normalized factor for each ROI. Q3 normalization was performed for all targets that were above the limit of quantification (LOQ) to reduce variance from segment size, segment cellularity and other technical factors.23

For our study, we used FFPE tissue samples from 56 patients with NSCLC who received immunotherapy. Four independent TMAs were constructed, each containing tumor tissue cores from these 56 patients. We made six independent cohorts as follows (table 1, figure 1). The CK-discovery cohort contained samples where the tumor tissue was intact and evaluable in all four TMAs. 136 patient cores (34 patients, four samples each) with adequate tumor consisted of CK-fourfold redundant cores. Then, the remaining 22 patient cores that were not present in all four TMAs created our CK-validation cohort. Similarly, CD45-fourfold redundant (22 patients, four cores each), CD45-validation cohort (31 patient samples) and CD68-fourfold redundant (12 patients, four cores each), CD68-validation cohort (43 patient cores) were created.

Two different analyses were performed. First, for the fourfold redundant TMA, we performed four independent Cox survival analysis and logistic regression (LR) models, respectively. Then we identified the genes that were significantly associated with outcome in each cohort and selected the ones that were common in all four cohorts. Finally, we performed Cox survival analysis and LR models in the validation cohort to see if the genes identified in the fourfold redundant TMA are significantly associated with the outcome in each compartment-specific validation cohort.

Statistical analysis

For quality assurance checking, normalization and batch effects adjustments, dimension reduction analysis using principal component analysis (PCA) was performed using The Dimension Reduction (V.1.1 updated March 2021) DSP DA script to control for batch variation. The SpatialDecon DSP (V.1.1 updated April 2021) script was used to estimate the abundance of mixed cell types within each ROI in this experiment. Linear mixed model was performed to identify differential gene expression between compartments. Pathway analysis (Reactome) was used to explore the unique landscape of gene expression across patients with DCB and compared with no durable benefit (NDB) post-ICI treatment. This was achieved by running Gene Set Enrichment Analysis (GSEA),24 using defined sets of genes per compartment statistically significantly associated with DCB. single sample GSEA scores are shown for the top key immunological pathways which exhibit diverse context-specific behaviors relative to DCB.

Continuous log-scaled data was used to identify significant associations between markers and treatment response in terms of OR using LR model. For the survival analysis, we first transformed the expression data to Z scores and calculated the HR using the Cox proportional hazards (Cox) survival model. Benjamini-Hochberg false discovery adjustment method was performed, considering the number of comparisons performed per compartment (tumor/CK, CD45, CD68). All hypothesis testing was performed at a two-sided significance level of α=0.05.

The association between the target expression and DCB was computed by LR using the R package “RegParallel” and “Biobase”. Survival analysis was performed using the R packages “RegParallel” and “Biobase”. To create the Kaplan-Meier plot representing the association of RPL13A and GNL3 expression with survival, samples were categorized as either “high” or “low” by using X-tile optimal cutpoint analysis.25 Reported unadjusted p values in the plot are associations of the model (log-rank test).

The software used for visualization of the data and statistical purposes was GraphPad Prism V.9.0 software for Windows (GraphPad Software, La Jolla, California, USA), and R studio V.1.4.1717. The optimal survival cutpoint analysis was performed using X-tile.26 In addition, the difference between GNL3 and RPL13A expression in NSCLC and normal tissue and the Cox-analysis between the two genes and non-ICI treated patients with NSCLC was studied by Tumor Immune Estimation Resource (TIMER V.2.0, http://timer.cistrome.org/).

Results

Our initial cohort included 56 NSCLC cases, each represented by four TMA cores in four master blocks. 54,000 transcript variables per ROI were generated, where 18,815 RNA markers, were measured in three molecularly defined tissue compartments (CK+tumor cells; CD45+leukocytes; CD68+macrophages). After filtering out RNA targets with an expression lower than the LOQ in at least 5% of segments, 8,580 RNAs from the 18,000 were evaluable for analysis (online supplemental table 3).

Biomarkers of resistance in the tumor compartment (CK+)

34 patients had adequate tumor in all four YTMAs tested. From YTMA1 121 RNAs out of 8,580 in the tumor (CK+) compartment were significantly associated with DCB (p<0.05). From YTMA2 217 RNAs; YTMA3 75 RNAs; YTMA4 205 RNAs. The intersection of the significant genes in the four independent TMAs containing tumor cores from the same 34 patients showed that increased expression of two RNA transcripts TFCP2L1 and CYBA was significantly associated with worse outcome (OR>1). Then we tested our validation cohort, which consisted of 22 patient cores that were not included in the fourfold redundant cohort. TFCP2L1 and CYBA were not associated with outcome in the validation cohort (online supplemental table 2).

We then sought to find potential markers that indicate resistance to immunotherapy regarding survival, irrespectively the tumor core was assessed. Taking the intersection of all the significantly associated RNAs with OS throughout the four YTMAs, RPL13A, GNL3, FAM83A, CYBA, ACSL4, SLC25A6, EPAS1, RPL5, APOL1, HSPD1, RPS4Y1 and ADI1 were correlated with shorter OS regardless the tumor core (HR>1, PLogRank<0.05). Two of them, ribosomal protein L13a (RPL13A) and GNL3 (both participating in ribosomal RNA (rRNA) processing in the nucleolus and cytosol) were also associated with shorter OS in the validation cohort (HR 2.48, PLogRank 0.02; HR 5.33, PLogRank 0.04) (figures 1 and 2figures 1–3, tables 2 and 3). Furthermore, RNA sequencing data from The Cancer Genome Atlas (TCGA) did not show a statistically significant correlation between RPL13A and GNL3 with survival in non-ICI treated patients with NSCLC but showed a higher expression in tumor compared with normal lung tissue (online supplemental figure 3 and table 3). Regarding PFS, CYBA and AP1M1 were associated with shorter PFS in all four TMAs of the 34 patients but not in the validation cohort.

Table 2

Cox proportional-hazards model analysis was performed for overall and progression-free survival. RPL13A and GNL3 were significantly associated with resistance to immunotherapy

Table 3

Overall survival for the 12 common genes (tissue microarrays 1–4) in the validation set; RPL13A and GNL3 are significantly associated with survival

Biomarkers of resistance in the immune cell compartments (CD45+ and CD68+ compartments)

22 patients had adequate lymphocytes (CD45+cells) in their tumor cores in all four YTMAs. From the intersections of the 94 (YTMA1), 93 (YTMA2), 79 (YTMA3) and 164 (YTMA4) significantly associated with DCB RNAs, NELFE was common in YTMA2 and YTMA3 and LIMS3 was common in YTMA3 and YTMA4. For the OS analysis, high secreted frizzled-related protein 2 (SFRP2, that modulates the influx from extracellular calcium in the B-cell receptor signaling pathway) in CD45 was significantly associated with shorter survival through all patients’ cores (PLogRank<0.05), but it did not validate in the independent cohort (PLogRank=0.54). Mitochondrial carrier 1 (MTCH1, that governs ferroptosis by triggering the FoxO1-GPX4 axis-mediated retrograde signaling) was associated with shorter PFS in YTMA2 and YTMA3 but not in the other YTMAs.

12 patients had adequate macrophages (CD68+cells) in their tumor cores in all four YTMAs. No common significantly associated with DCB RNAs were found across all four YTMAs of the 12 patients in CD68 compartment. When we ran the survival analysis, Rho GTPase-activating protein (ARHGAP, a negative regulatory factor of Rho family proteins27) in the CD68 compartment was significantly associated with shorter OS in YTMA 1, 2 and 3 and in the validation cohort (HR 10.35, PLogRank 0.04). PNN interacting serine and arginine rich protein in the CD68 compartment was significantly associated with favorable PFS in YTMA 1, 2 and 3 and in the validation cohort (HR 0.23, PLogRank 0.03).

Differential expression analysis and pathway analysis

Differential expression analysis among compartments showed differences between tumor and CD45, CD68 compartments. 1,658 and 3,083 out of 8,580 genes were significantly differentially expressed between tumor (CK), CD45 and CD68, respectively (P adj<0.05). RPL13A and GNL3 expression was higher in the tumor (CK) compartment compared with the immune compartments (online supplemental figure 4).

The panel of genes profiled has a strong focus around capturing biological signaling along canonical signaling pathways (Reactome pathway database) and cell-intrinsic signaling from tumor and immune cells. We tried to leverage foundational pathway interrogation tools with the data to capture spatially resolved pathway interactions and signaling within the FFPE tissues of ICI treated lung cancer patients. In CK compartment interferon (IFN) signaling (Normalized Enrichment Score, NES=3.84, adjusted p (adj p)=0.0030), cell cycle checkpoints (NES=3.36, adj p=0.0030), antigen processing-cross presentation (NES=3.16, adj p=0.0030) were among the top enriched pathways associated with DCB to immunotherapy (online supplemental table S4). Interestingly metabolism of steroids pathways were enriched for resistance to immunotherapy in all three compartments: CK (NES=−2.46, adj p=0.0058), CD45 (NES=−2.4, adj p=0.0045) and CD68 (NES=−2.18), adj p=0.0058) (online supplemental figure 5). Steroid regulation of beta-adrenergic receptor might have a links to T-cell exhaustion and response against cancer.28 Ribosome protein pathway is a regulator of lipid homeostasis and metabolism of steroids and might influence the efficacy of immunotherapy.29 While these pathways are enriched in non-responders, future efforts will determine if gene subsets might have value as compartment specific tumor signatures for outcome.

Discussion

In this study, we used DSP technology as a high-plex screening tool to identify biomarkers of resistance to ICI treatment. We leveraged the spatial information separating tumor from host response to identify protein-coding genes that are associated with outcome. To address the intratumor heterogeneity we analyzed each sampling site’s genes independently and then took the intersection of the ones that were identified as significantly associated with outcome to increase the reproducibility of the observation and robustness to spatial heterogeneity. We noticed that in the tumor (CK+) compartment 12 protein-coding genes (RPL13A, GNL3, FAM83A, CYBA, ACSL4, SLC25A6, EPAS1, RPL5, APOL1, HSPD1, RPS4Y1 and ADI1) were correlated with survival throughout the different tumor cores. In the independent validation cohort high expression of two, ribosomal protein L13a (RPL13A) and G protein nucleolar 3 (GNL3), were also associated with worse post-ICI treatment outcome. These RNA-binding proteins coding genes were correlated with shorter survival in two independent ICI-treated NSCLC cohorts. We did not find any association of RPL13A and GNL3 with survival in TCGA non-ICI treated cohort. Thus, RPL13A and GNL3 is indicative but not prognostic for survival.

rRNA processing is a highly active and essential process in lung cancer cells.30 Ribosome biogenesis has attracted increasing attention in molecular research over the past decade due to its beneficial role in controlling the translation of all proteins in the cell, and thus governing cell growth and proliferation.31 The immunological value of the regulators of rRNA modification is under investigation in different cancer types.32 33 As an example, increased expression of RPL14 and RPL28 induced by MYC overexpression in vivo is a major event in lymphoma progression, leading to impaired translation of the cell cycle regulator CDK11, genomic instability, and defective mitosis that directly contribute to lymphomagenesis.34 35 Increasing evidence has demonstrated that RNA modifications have indispensable roles in inflammation, innate immunity, and antitumor activities, among which the most extensive and in-depth research has been on m6A. Zhang and colleagues demonstrated that modification of m6A methylation is involved in the regulation of the microenvironment of gastric cancer, and has a guiding role in immunotherapy.36 In colon cancer a tumor immune infiltration-associated small nucleolar RNA model (which included GNL3) has been proposed as a potential biomarker of immunotherapy response.37 In NSCLC, RNA-binding motif 10, a tumor suppressor protein, can directly target c-Myc for degradation and reduce its cancer-causing effects by binding with RPL5 and RPL11.38 Inhibition of rRNA transcription and synthesis has been shown to induce apoptosis and cell cycle arrest as key modulators of both the carcinogenic process and the immune response30 39 and, thus, it seems to influence the efficacy of immunotherapy.

Preclinical data in ovarian cancer showed that antitumor activity of OTX015A, a novel and selective BET inhibitor,40 was linked to a marked downregulation of GNL3 gene.41 GNL3, an evolutionarily conserved stem cell gene is influencing cell proliferation.42 Nucleostemin (NS; a product of the GNL3 gene) was shown to be involved in DNA repair pathways and modulate PD‐L1 and cytokine expression in different types of cancer.43 44 In addition, the overexpression of RPL13 promoted the induction and activation of the promoters of the nuclear factor-κB and IFN-β genes, and the expression and protein secretion of the antiviral factor IFN-β and proinflammatory cytokine interleukin-6.45 RPL13A can also act as a negative regulator of the inflammatory response.46 In the late stage of IFN-γ stimulation, the IFN-γ-activated inhibitor of translation complex formed by RPL13A, glyceraldehyde-3-phosphate dehydrogenase, glutamyl-prolyl-tRNA synthetase, and NS1-associated protein 1, inhibits the translation of related inflammatory genes.

One of the limitations of this study is that we analyzed tumor samples collected from retrospective immunotherapy-treated cohorts from a single institution, not from randomized clinical trials. Hence, data collection occurred before the study’s initiation, and so we did not have control over the sample size or the allocation of participants, which are key components for conducting a power calculation. Larger cohorts of patients with NSCLC treated with ICI are warranted to further investigate the role of GNL3 and RPL13A expression on ICI resistance. For this study FFPE tissue samples were used. Work from our laboratory has shown that the rate of antigenicity loss is biomarker specific and should be considered when archived tissue is used.47 48 41 out of 56 biopsies used for this study are dated after 2016 so the time interval from collection to study initiation (2022) was less than a decade, that is generally considered as acceptable (online supplemental table S5).49 Another limitation is that we did not fully validate an identical cutpoint in the training and validation sets for optimal patient stratification due to the lack of true standardization of assessment by DSP. Future studies addressing the predictive value of rRNAs (RPL13A and GNL3) will require the validation of an optimal and reproducible cutpoint in larger well-powered cohorts. Although there is substantial data to support the role of these rRNAs as negative regulators of the immune response, further functional studies are needed to understand the interplay between GNL3 and RPL13A expression levels and mechanisms of immune evasion.

In conclusion, diverse mechanisms of immunotherapy resistance have been characterized, and more continue to be uncovered.50 51 rRNAs have multiple functions as systemic regulators in biological processes and may also play a crucial role in immunotherapy resistance.52 Future explorations are warranted to discover more immune-related rRNAs and elucidate their common mechanisms of immune regulation NSCLC, thus giving more direct indications to deal with immune resistance.

Data availability statement

Data are available upon reasonable request.

Ethics statements

Patient consent for publication

Acknowledgments

The authors thank Lori A Charette and the staff of Yale Pathology tissue services for expert histology services, G Wang, J Heltke, C Castaldi and the stuff of Yale Center for Genome Analysis (YCGA) for sequencing clinical samples services. DLR, is supported by the Yale SPORE in Lung Cancer and the Yale Cancer Center CCSG and DLR was supported by a sponsored research agreement from Konica Minolta for part of this effort. MKM received a scholarship from the Hellenic Society of Medical Oncologists (HESMO).

References

Supplementary materials

Footnotes

  • X @mkmoutafi, @tznaung, @rolando_milian, @IVathiotis, @AthAngelakis

  • Contributors Conceptualization: DLR, MKM. Methodology: MKM, DLR, TNA. Software: RGM, AA, LS. Validation: TNA, MKM. Formal analysis: MKM, RGM, AA, LS. Investigation: MKM, DLR, KMB, VX. Data curation: MKM, IAV, NG, KAS. Writing (original draft): MKM, DLR. Writing (review and editing): MKM, DLR, KMB, VX, LS, KAS. Visualization: MKM. Supervision: DLR. Project administration: DLR. Funding acquisition: DLR. Content guarantor: DLR.

  • Funding This work was supported in part by Yale SPORE in Lung Cancer, P50 CA196530

  • Competing interests DLR has served as an advisor for AstraZeneca, Agendia, Amgen, BMS, Cell Signaling Technology, Cepheid, Danaher, Daiichi Sankyo, Genoptix/Novartis, GSK, Konica Minolta, Merck, NanoString, PAIGE.AI, Roche, and Sanofi. Amgen, Cepheid, NavigateBP, NextCure, and Konica Minolta fund or have funded research in DLR’s laboratory. KAS has served as a Consultant or advisor: AbbVie Inc.; Agenus Inc.; AstraZeneca; Bristol Myers Squibb; CDR-Life Inc.; Clinica Alemana de Santiago; EMD Serono Inc.; F.Hoffmann-La Roche; Genmab A/S; Indaptus Therapeutics; Janssen Pharmaceuticals, Inc.; Merck Sharpe & Dohme; Moderna, Inc.; Molecular Templates, Inc.; OnCusp Therapeutics; Parthenon/Incendia Therapeutics; Repertoire Therapeutics; Sanofi; Sensei Biotherapeutics; Shattuck Labs, Inc.; and Takeda Pharmaceutical Company Limited.

  • 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.