Single-cell-eQTL mapping in circulating immune cells reveals genetic regulation of response-associated networks in lung cancer immunotherapy

While immune checkpoint inhibitors (ICIs) are adopted as standard therapy for advanced non-small cell lung cancer (NSCLC), genetic determinants of response heterogeneity remain elusive [1]. As most hematopoietic lineages undergo dynamic changes during tumor pathogenesis and immunotherapy [2], elucidating how germline variants modulate their transcriptomes is critical. Expression quantitative trait loci (eQTL) analysis, especially integrated with single-cell RNA sequencing (scRNA-seq), enables gene regulation mapping at single-cell resolution [3, 4]. Detailed methodologies are described in the Supplementary Materials.

To investigate how germline variants shape immune gene regulation during ICI treatment, we performed single-cell-eQTL (sc-eQTL) analysis and transcriptomic network profiling. Peripheral blood mononuclear cells (PBMCs) were collected from 73 NSCLC patients treated with anti-programmed cell death protein-1 (PD-1) or programmed death-ligand 1 (PD-L1) therapy, at both baseline and 1-5 weeks post-treatment (Figure 1A, Supplementary Table S1). By integrating scRNA-seq with SNP array data, we analyzed cell-type-resolved sc-eQTLs and gene networks (Figure 1A-B).

After quality control and pseudobulk aggregation, we identified 9,147 eQTL pairs—expression-regulating SNPs (eSNPs) linked to 3,616 blood expression-regulated genes (eGenes)—across eight immune cell clusters and treatment conditions (Figure 1B-C, Supplementary Figure S1A, Supplementary Table S2). Consistent with previous studies [3, 4], eGene counts correlated with cell abundance, and eSNPs were enriched in regulatory elements (Supplementary Figure S1B-D). Multiadaptive shrinkage [5] revealed distinct cell-type- and treatment-dependent regulation, including 245 treatment-specific eQTLs (Supplementary Figure S2A-D, Supplementary Tables S3-S4). For instance, tumor necrosis factor (TNF) was regulated in monocytes post-treatment (posterior β = 1.17), while TNF receptor 1A (TNFRSF1A) was baseline-regulated in CD8⁺ T cells (β = 1.10), indicating genetic variants may shape immune gene expression during ICI therapy (Figure 1D). Additional examples include key cytotoxic mediators perforin 1 (PRF1) and granzyme B (GZMB) in baseline CD8⁺ T cells (Figure 1D, Supplementary Figure S2E).

To validate our findings, we conducted two complementary analyses. First, colocalization analyses with genome-wide association study (GWAS) loci for autoimmune and blood traits showed overlaps (PP.H4 > 0.6), suggesting possible shared regulatory mechanisms (Supplementary Figure S3, Supplementary Tables S5-S6). Second, comparison with external eQTL studies showed our study-specific eQTLs, hereafter referred to as lung cancer-specific eQTLs, were enriched in cancer- and immune response-related pathways (Figure 1E, Supplementary Figure S4A), reflecting chronic immune adaptation. Nevertheless, most regulatory patterns seen in healthy donors from the 1M-scBloodNL study [3] were preserved (Supplementary Figure S4B, Supplementary Tables S7-S8).

Beyond eQTL mapping, we used weighted gene co-expression network analysis (WGCNA) to build co-expression networks for each major immune cell type [6], enabling comparison of cell-type-specific regulation and transcriptomes (Supplementary Tables S9-S10). Among the identified modules, the CD8-brown module showed elevated baseline activity in the non-durable clinical benefits (NCB) group (Figure 1F, Supplementary Figure S5A). This module was cytotoxic-gene-rich (PRF1, apolipoprotein B mRNA-editing enzyme, catalytic subunit 3G [APOBEC3G], and GZMB) and highly expressed in the NCB group, particularly in differentiated CD8⁺ T subclusters (Supplementary Figure S5B-D). Its activity was supported by external blood datasets from healthy individuals, tumor patients, and ICI-treated tumor datasets (Supplementary Figure S5E-F).

To explore potential regulators, we applied single-cell regulatory network inference and clustering (SCENIC) [7], identifying eomesodermin (EOMES) and t-box transcription factor 21 (TBX21) as putative regulators linked to CD8⁺ T cell differentiation and exhaustion [8]. Their regulon activity aligned with both the CD8-brown module (Figure 1G, Supplementary Figure S6A-B), and the abundance of effector memory CD8⁺ T (CD8⁺ TEM) cells (Supplementary Figure S6C-D). Core genes, including EOMES, TBX21, and interleukin-2 receptor, beta subunit (IL2RB), were up-regulated in NCB compared to durable clinical benefit (DCB) group (Supplementary Figure S6E). Notably, TBX21 and EOMES were not eGenes in our analysis.

Building on these findings, we employed a graph neural network (GNN)-based framework to refine the CD8-brown module by integrating protein-protein interaction (PPI) networks, gene ontology (GO) annotations, and co-expression data (Supplementary Figure S7A-B). The refined subnetwork captured a more coherent CD8⁺ T cell differentiation program and showed stronger topological connectivity than the original WGCNA-derived gene set in Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; node density: 0.333 vs 0.183; Supplementary Figure S7C-D, Supplementary Table S11) [9]. Pathway enrichment for canonical pathways in CD8⁺ T cells supported this refinement (Supplementary Figure S7E-F).

This module was well-conserved across external scRNA-seq datasets, including the tumor microenvironment (TME) of ICI-treated NSCLC patients (Supplementary Figure S8) [10] and blood from cancer patients (Supplementary Figure S9). Its activation was significantly enhanced in non-responders, particularly in differentiated CD8⁺ T cell subclusters like CD8⁺ TEM (Figure 1H, Supplementary Figures S8-S9, Supplementary Table S12), suggesting that this systemic cytotoxic signature is recapitulated in the tumor microenvironment and may reflect ICI response-associated immune status.

Given its strong association with ICI response, we evaluated whether the refined core genes from the CD8-brown module could explain survival outcomes. We derived a core gene score from 15 GNN-prioritized CD8-brown genes with high network centrality (Supplementary Methods, Supplementary Figure S7, Supplementary Table S13), which strongly correlated with module activity (Supplementary Figure S10A). Using grid search and bootstrapping (Supplementary Figure S10B-C), we determined a threshold to classify patients: those with >75% of CD8⁺ T cells with high CD8-brown core gene score (>0.36) displayed significantly shorter overall survival (OS, P = 0.035) and progression-free survival (PFS, P = 0.018; Figure 1I, Supplementary Figure S10B-C). The survival association remained significant in a multivariate Cox model adjusting for clinical covariates (Supplementary Figure S10D). Overall, our analysis suggests that the TBX21-EOMES regulatory axis may drive the CD8-brown module, representing CD8⁺ T differentiation with cytotoxic eGenes like PRF1 and GZMB, which in turn stratify poor responders to ICI.

Finally, to evaluate how germline variants influence the immune network associated with ICI response, we examined the distribution of eQTLs within co-expression and regulatory networks (Figure 1J-L, Supplementary Figure S11A-C). eGenes were enriched in central network positions (Supplementary Figure S11B); however, no single eQTL–including PRF1, GZMB or TNFRSF1A–exerted a dominant regulatory effect on the CD8-brown module. Instead, genes with higher centrality or WGCNA-assignment were linked to smaller eQTL effect sizes (Figure 1K-L), suggesting functional constraints on key immune differentiation pathways.

Despite these insights, several limitations exist. First, dichotomizing responses into DCB and NCB may obscure heterogeneous outcomes, while the treatment-induced transcriptomic changes appeared smaller than inter-sample heterogeneity (Supplementary Figure S12). Second, although monocyte-specific eQTL were detected and monocytes-ICI interactions are well-known, we excluded them from downstream analysis due to minimal contribution to ICI-associated co-expression modules (Supplementary Figure S13). Third, our focus on cell-type- and condition-specific eQTLs and co-expression networks may overlook intercellular interactions and trans effects influencing immune function. Lastly, although the CD8-brown core gene score stratified survival in our cohort (c-index = 0.68-0.71 for PFS/OS, Supplementary Figure S10D), its predictive value requires validation in independent datasets.

Taken together, our study identified 3,616 blood- and 702 lung cancer-specific eGenes via sc-eQTL mapping and revealed a transcriptomic module (CD8-brown) associated with CD8⁺ T cell differentiation and non-responder for ICI response. The limited influence of eQTLs on co-expression networks suggests functional constraints in the immune transcriptome. By linking genetic variation to cytotoxic network activity and clinical outcomes, our analysis provides a framework to understand systemic immunity under ICI treatment in metastatic NSCLC. Future research involving larger cohorts and experimental validation of key regulatory variants—such as those affecting PRF1 and GZMB—could clarify causal mechanisms and refine personalized ICI strategies.

Hyungtai Sim: Data curation; formal analysis; validation; investigation; writing—original draft. Geun-Ho Park: Data curation; formal analysis; investigation. Woong-Yang Park: Investigation. Se-Hoon Lee: Conceptualization; investigation; writing—original draft. Murim Choi: Conceptualization; investigation; writing—original draft.

The authors declare no competing interests.

This work was supported in part by grants from the Korean Research Foundation (NRF-2021R1A2C3014067, NRF-RS-2023-00207857 to Murim Choi, and NRF-2020R1A2C3006535, NRF-RS-2025-00519956 to Se-Hoon Lee), the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (HR20C0025 to Se-Hoon Lee), and Future Medicine 20*30 Project of the Samsung Medical Center (SMX1230041 to Se-Hoon Lee).

This study investigated human samples and was given permission by the Samsung Medical Center Institutional Review Board (No. 2018-04-048, 2022-01-094). All participants underwent an informed consent process before being enrolled in the study.