Integrated Single-Cell Network and LLM-Based Analysis Reveals Activation and Recognition Mechanisms in γδ T Cells in Human Colorectal Cancer

Introduction: : γδ T cells, a unique T cell subset enriched in the colon epithelium, patrol the mucosal barrier for signs of infection and malignancy and contribute to anti-tumor immunity in colorectal cancer (CRC). They are promising candidates for immunotherapeutic strategies because they are naturally present in the tumor microenvironment, exhibit strong tissue infiltration, and recognize targets in an MHC-independent manner, enabling potential allogeneic transfer. However, despite this promise, the efficacy of γδ T cells in tumor elimination can still be improved: our previous work integrating single-cell RNA sequencing (scRNA-seq) data from multiple CRC studies revealed diverse γδ T cell states, including a large population of unexhausted, tissue-resident bystanders expressing ITGAE and ITGA1, alongside highly cytotoxic effectors producing IFNG. Unfortunately, the mechanisms underlying the activation and tumor recognition of these bystanders remain poorly understood. Identifying the molecular cues that drive their activation could support early tumor elimination and prevent recurrence after surgery. Here, we combined multi-study single-cell transcriptomic integration, gene regulatory network inference, and large language model (LLM)-based in silico gene perturbation to uncover the molecular drivers of γδ T cell activation in CRC.

Materials and

Methods: : We hypothesized that γδ T cells are activated by ligands present in the tumor microenvironment that engage their surface receptors. To identify candidate ligands driving the shift from tissue-resident memory (TRM)-like bystanders to highly cytotoxic effectors, we first applied SCENIC to infer a γδ T cell–specific gene regulatory network (GRN), addressing the limitation that existing signaling networks are not tailored to this cell type. We then integrated the SCENIC-derived GRN into NicheNet, a computational framework that performs a random walk on the known intracellular signaling network to infer ligands responsible for observed differential gene expression. Ligand predictions (ligands with AUROC>0.5 were kept) were based on differentially expressed genes (Log2FC > 1, adjusted p < 0.05) between TRM-like and effector γδ T cell subsets.

To identify intracellular genes capable of causally driving this state transition, we fine-tuned Geneformer, a large language model (LLM) pretrained on 95 million human single-cell profiles that enables in silico perturbation analysis by modeling how gene expression changes shift a cell’s latent embedding in relation to different cell states. We fine-tuned Geneformer using our tumor-infiltrating γδ T cell data, freezing six layers and applying 100 iterations of hyperparameter optimization. We then performed in silico overexpression perturbations across 5,000 less-activated γδ T cells and retained genes whose perturbation shifted cells toward effector-like embeddings (FDR < 0.05, cosine shift > 0.003, and expression in ≥50 perturbed cells). Finally, we re-analyzed these Geneformer-predicted genes using the γδ T cell–refined NicheNet model to cross-validate the predicted upstream ligands.

Results, Conclusions, and Discussions:: SCENIC analysis predicted over 13,200 novel transcription factor (TF)–target gene interactions in CRC-infiltrating γδ T cells. Integration of these interactions into the known signaling network did not compromise the model’s ability to predict perturbed ligands in reference data. Our analysis predicted that ligands such as PD-L1 and LGALS9 contribute to γδ T cell exhaustion in CRC, while IL4, IL1B, TNFSF9, and IL-15 were predicted to enhance γδ T cell activation—though this activation remains influenced by suppressive ligands including TGFB1 and PD-L1.

Using the annotated CRC whole-tissue scRNA-seq reference dataset from Joanito et al. (2022), we further assessed the enrichment of top γδ T cell phenotype-shifting ligands across various cell types within the CRC microenvironment. Our findings indicate that in both MSS and MSI-H CRC, monocytes/classical dendritic cells (McDCs) have the highest potential to promote γδ T cells' transition from a less activated, TRM-like state to a more activated, effector-like state.

In silico overexpression perturbations using the fine-tuned Geneformer model identified several genes—CCR9, NCR2, KLRK1 (NKG2D), KLRC3 (NKG2E), NR4A3, and MMP9—whose activation causally shifted tissue-resident bystander γδ T cells toward an effector-like state. Importantly, most of these genes were not differentially expressed between effector and TRM-like cells, indicating that Geneformer identified novel regulators beyond conventional differential expression analysis.

The γδ T cell–refined NicheNet was then re-applied to the Geneformer-predicted genes whose overexpression promoted an effector-like phenotype. The predicted upstream ligands of these genes included IL15 and TNFSF9, which were also among the top ligands predicted to explain the differentially expressed gene signatures.

Together, these findings highlight IL15 and TNFSF9 as promising candidates for promoting γδ T cell activation in CRC. In addition, NCR2 and NKG2E emerged as potential receptors used by γδ T cells to recognize CRC cells. These insights inform the design of engineered γδ T cell therapies to enhance anti-tumor function and support early tumor elimination and post-surgical surveillance.

Acknowledgements and/or References (Optional)::