Elucidating the Mechanisms of SA–4–1BBL-Mediated Cancer Immunoprevention Through Advanced Informatics Approaches

TL;DR

This study uses bioinformatics to uncover how SA–4–1BBL prevents cancer by triggering a balanced immune response.

Contribution

The paper reveals the distinct immune mechanisms of SA–4–1BBL compared to antibodies in cancer immunoprevention.

Findings

SA–4–1BBL activates early effector immune cells at the injection site and sustained regulation in lymph nodes.

SA–4–1BBL induces selective adaptive immunity and modulates the CD151–TGF-β axis for balanced immune response.

3H3 antibody activates broader innate inflammatory pathways, unlike the more targeted SA–4–1BBL response.

Abstract

Cancer immunoprevention leverages the immune system’s surveillance mechanisms to mitigate tumor development. Vaccines that constitute a tumor antigen and an immune adjuvant are perceived as immunoprevention modalities. However, relevant tumor antigens are unknown for non-viral cancers, which constitute most human cancers. Our group has recently shown that SA–4–1BBL, a novel agonist of CD137 receptor, but not antibodies, shows immunoprevention efficacy against various tumors. Advanced bioinformatics analyses of bulk RNA-seq data were conducted to elucidate mechanisms underlying cancer immunoprevention. Mice received subcutaneous injections of SA–4–1BBL or agonistic 3H3 antibody, and the injection-site tissue (IS) and draining lymph nodes (LN) were analyzed for differential gene expression. SA–4–1BBL induced a compartmentalized and temporally dynamic immune program characterized by early effector activation at IS and sustained immune regulation in draining LN. K-means clustering of 4564 DEGs identified eight functionally distinct clusters. IS-enriched clusters contained activation genes for CD4+ T and NK cells, including Cd28, Klra1, Cd4, Cd40, and Cd40l, while LN clusters were enriched for regulatory genes (Tnfaip3, Irf5, Col1a2) that ensure immune priming and homeostatic restraint for a balanced response. SA–4–1BBL generated a more selective and durable activation of adaptive immunity, TCR signaling, Th1/Th2 differentiation, and NK cytotoxicity. 3H3 activated broader innate inflammatory programs, including Toll-like receptor and neurodegeneration-linked pathways. IMPRes analysis showed that SA–4–1BBL activates sequential immune-regulatory circuits centered on Stat1, Cd247, and Ifng and modulates the CD151–TGF-β axis. These findings demonstrate that SA–4–1BBL elicits a balanced immune response, ensuring both safety and efficacy in preventing cancer development.

Figures

Full text

1. Introduction

The immune checkpoint stimulator 4–1BB (CD137; TNFRSF9), as a type I transmembrane protein, is one of the most potent members of the TNF receptor superfamily. It shows induced expression in activated CD8^+^ and CD4^+^ T, NK, and NKT cells, and constitutive expression in dendritic cells, neutrophils, monocytes, and regulatory T cells [1,2,3,4,5]. 4–1BB signaling enhances the proliferation and survival of T cells, the acquisition of effector function, and the establishment of long-term memory. Although 4–1BB signaling primarily targets CD8^+^ T cells, it also enhances the cytotoxic function and cytokine production of NK cells, particularly increasing IFN-γ, which contributes to their anti-tumor response [6,7,8]. Signaling via 4–1BB also reverses T cell exhaustion and overcomes T regulatory cell suppression, making this pathway an attractive target for cancer immunotherapy [2,9,10]. Agonistic antibodies (Ab) to 4–1BB have shown anti-tumor efficacy in preclinical models, leading to their evaluation in clinical trials. However, both preclinical and clinical studies have demonstrated that agonistic Abs can cause significant toxicity manifested by cytokine storm, alteration of immune cell trafficking, and depletion of various immune cells, including B, NK, and CD4^+^ T cells [11,12,13]. The first clinical trial with an agonistic monoclonal Ab, urelumab (BMS-663513), developed by Bristol–Myers Squibb for immunotherapy (NCT00309023) against melanoma and renal cell carcinoma, showed significant liver toxicity [14]. Another Ab, utomilumab (PF-05082566), developed by Pfizer in combination with an anti-PD-1 antibody, pembrolizumab, although it showed reduced toxicity, had moderate therapeutic efficacy in patients with advanced solid tumors [15]. Lack of efficacy and observed toxicity may not be a bona fide feature of 4–1BB signaling but rather an attribute of Abs that cannot recapitulate the signal delivered by natural ligands [16,17].

Effective 4–1BB signaling requires cross-linking by its membrane-bound dimeric (mouse) or trimeric (human) ligand (4–1BBL) expressed on antigen-presenting cells, including monocytes, macrophages, and dendritic cells [18,19,20,21,22]. Soluble 4–1BBL lacks immune-stimulating function as it cannot crosslink 4–1BB receptor on immune cells, limiting its use for cancer immunotherapy [23,24,25]. To overcome this challenge, we previously reported the generation of a novel recombinant form, SA–4–1BBL, containing the extracellular domain of mouse ligand chimeric with a modified form of core streptavidin [13]. This novel agonist forms tetramers and oligomers owing to the structural features of streptavidin and has robust immunoregulatory functions [26,27]. Importantly, SA–4–1BBL does not elicit systemic proinflammatory cytokine storm and severe liver toxicity as reported for agonistic 4–1BB Abs [13]. As an adjuvant component of tumor-associated antigens-cancer vaccines, SA–4–1BBL showed robust therapeutic efficacy in various preclinical tumor models without reported toxicity [13,28,29,30,31,32]. Surprisingly, SA–4–1BBL as a single agent demonstrated prophylactic efficacy against transplantable tumors in mice. In marked contrast, agonistic 4–1BB Abs lack such efficacy [33]. The primary objective of this study is to investigate the mechanisms underlying the cancer immunoprevention efficacy of SA–4–1BBL using bulk RNA sequencing (RNA-seq) and advanced translational bioinformatics approaches.

Bulk RNA-seq enables comprehensive analysis of differential gene expression (DGE) across various tissues to assess treatment-specific effects [34]. For example, analysis of DEGs at the injection site draining lymph nodes (LN), which serve as the nucleus for immune regulation, can identify pathways important for immune surveillance of cancer [35]. Such analyses can reveal how SA–4–1BBL modulates various immune cells, including T cells, NK cells, and myeloid populations, which are critical for its broad immunomodulatory effects [36]. Similarly, analyzing DEGs at the injection site (IS), a primary locale of antigen presentation and innate immune activation, can uncover early signaling events that shape downstream immune effector responses. Insights from spatial transcriptomics can elucidate how SA–4–1BBL coordinates immune programming across different tissues to enhance tumor immunosurveillance, which will facilitate the clinical development of this agent for cancer immunoprevention [37].

2. Materials and Methods

2.1. Mice and Reagents

C57BL/6 mice (6–10 weeks) were purchased from Jackson Laboratory (Jackson, ME, USA) and maintained and treated according to protocols approved by the Institutional Animal Care and Use Committee. SA–4–1BBL recombinant protein was produced and characterized for structure and function in our laboratory per published protocols [14,34,38]. An agonistic Ab (3H3) to CD137 was purchased from BioXcell BE0239, Clone 3H3, Lebanon, NH, USA.

2.2. Study Design

C57BL/6 mice were treated subcutaneously with 100 mg of SA–4–1BBL or equimolar of an agonist Ab (3H3) either once or twice, two weeks apart. The dose of SA–4–1BBL was based on prior studies demonstrating its robust cancer immunoprevention efficacy [36]. Mice were euthanized 3 or 7 days post last injection of SA–4–1BBL to capture changes in innate and adaptive immune pathways, respectively. Tissues from the injection site (skin) and draining lymph nodes were collected for analysis. Saline-treated mice served as controls. Each condition included four biological replicates, for a total of 32 mice (Figure 1). Lymph nodes and skin tissues were rinsed with PBS and stored at −80 °C in TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) to extract total RNA.

2.3. RNA-Seq

Tissues in TRIzol were homogenized using a handheld microtissue homogenizer (VWR200 Homogenizer, VWR International, LLC, Radnor, PA, USA), kept at room temperature for 5 min. After adding chloroform, the homogenate was vortexed for 30 s, incubated at room temperature for 5 min, and spun at 12,000× g for 15 min at 4 °C. The aqueous layer was transferred into a new sterile RNase-free 1.5 mL Eppendorf tube, followed by adding an equal volume of 100% ethanol to proceed with total RNA extraction. The RNeasy Mini Kit (Qiagen, Hilden, Germany) was used to isolate high-quality total RNA per the manufacturer’s instructions. RNA was stored in RNase-DNase-free water (Invitrogen, Carlsbad, CA, USA) and assessed for quantity and purity using Nanodrop One (Thermo Fisher Scientific), Agilent Bioanalyzer (Agilent, Santa Clara, CA, USA). Batches of RNA were shipped to LC Sciences (Houston, TX, USA) for poly(A) RNA sequencing. A cDNA library was prepared, followed by Illumina TrueSeq-standard sample preparation protocol, and paired sequencing was performed on the Illumina NovaSeq 6000 (San Diego, CA, USA) sequencing system. Across all samples, the average read count per sample ranges from ~17 million to ~25 million reads.

2.4. Bulk RNA-Seq Quality Control and Differential Gene Expression Analysis

Quality assessment of raw RNA-seq reads was performed using FastQC (v0.11.9) [39], with results consolidated via MultiQC (v1.11) [40]. Reads were then trimmed using Trim Galore (v0.6.7) with cutadapt [41] to remove Illumina adapters, ambiguous nucleotides, and sequences < 20 bp, applying a Phred quality score of 20 for 99% sequencing accuracy. Trimmed reads were aligned to the mouse genome (GRCm39.107) using HISAT2 (v2.2.1) [42], achieving a 97% alignment rate (Supplementary Table S1). Aligned files were converted to sorted BAM format with Samtools (v1.14) [43], and gene expression levels (FPKM values) were quantified with Cufflinks (v2.2.1). Differential gene expression analysis was performed using Cuffdiff (V2.2.1), which estimates transcript- and gene-level expression variability using a beta negative binomial modeling framework with replicate-aware variance estimation. This approach does not implement linear modeling or empirical Bayes variance moderation. These FPKM values enabled PCA plot construction for evaluating intrinsic sample differences. Prior to PCA, FPKM expression values were log2-transformed and scaled to reduce the mean–variance dependency of RNA-seq data and to ensure equal contribution of genes across samples in the principal component space.

In our study, the Cuffdiff tool was utilized to compute gene expression differences across 24 selected comparisons out of 49 total comparisons from 18 samples with a total of 62 replicates (Supplementary Table S2). DEGs were identified based on a minimum two-fold change (absolute log_2_ fold change ≥1) and a statistically significant q-value (q-value ≤ 0.05) in each selected comparison, extracted from the Cuffdiff results. Here, a “sample” refers to a unique mouse × tissue (IS or LN) × endpoint (D3/D7) × injection regimen (1X/2X) × treatment combination. Biological replicates are independent mice within the same condition. Technical replicates (if applicable, e.g., repeat library prep or resequencing) were handled by merging reads as replicates and are explicitly labeled in Supplementary Table S2. The phrase “18 samples” refers to 18 RNA-seq libraries, and “63 replicates” refers to the total libraries, including technical replicates across all tissues/conditions. From these libraries, we evaluated 49 possible contrasts and selected 24 primary contrasts for downstream integrative analyses (Supplementary Table S2) to focus on treatment, tissue, time, and regimen effects while maintaining interpretability.

2.5. Candidate Prioritization and Hypothesis Generation

To further prioritize our DEGs list for drawing inferences, we compared the absolute difference in Log_2_ Fold Change (Log_2_FC) between the “SA–4–1BBL vs. Naïve” and “3H3 vs. Naïve” comparisons, by defining a “trend parameter”, to keep track of whether the DEGs from the two comparisons were either both up-regulated (both up), both down-regulated (both down), or had one up-regulated and other down-regulated (opposite). To lie under any one trend criterion, the “SA–4–1BBL vs. Naïve” and “3H3 vs. Naïve” comparisons should show an absolute difference in Log_2_ Fold Change (Log_2_FC) ≥ 1 in both comparisons.

2.6. Distance Matrix and Barcode Strategy

The Euclidean Distance Matrix (EDM) was used independently to capture pairwise distances between samples based on DEG profiles, revealing clustering patterns across treatment, tissue, and time. Separately, a treatment-specific barcode strategy was employed to classify DEGs based on their expression trends in response to 3H3 or SA–4–1BBL treatments, enabling fine-grained annotation across multiple experimental conditions. This 6-digit barcode classification was derived from a total of 24 comparisons, eight each for 3H3 vs. Naïve, SA–4–1BBL vs. Naïve, and SA–4–1BBL vs. 3H3 as represented in the Treatment column of Figure 1C, where each comparison captures distinct DEG contrasts across treatment, tissue type, time point, and injection number. Categorized DEGs using a binary code (0 and 1), with representation as 0–1 (more prominent in SA–4–1BBL), 1–0 (more prominent in 3H3), or 1–1 (equally expressed). Treatment specificity was determined based on absolute Log_2_ fold change (Log_2_FC) values.

2.7. K-Mean Clustering and Random Forest Model-Based Ranking

K-means clustering was applied to systematically stratify differential gene expression patterns across experimental conditions. Clustering was performed on differentially expressed genes (DEGs) derived from 24 pairwise comparisons using the Scikit-Learn implementation in Python (V3.10) (Supplementary Figures S3 and S4). Genes were grouped based on their normalized expression profiles across treatments, tissues, and time points to identify coherent temporal and regimen-dependent expression trends.

The optimal number of clusters (K) was determined using the elbow method, in which the total within-cluster sum of squares (WCSS) was calculated across a range of K values (typically 1–10) to identify the point of diminishing returns. For validation and consistency, WCSS calculations were additionally assessed using the K-means function in R.

In parallel, a random forest (RF) model was employed as a supervised, non-linear feature-selection approach to rank genes based on their relative importance in distinguishing experimental conditions. The RF analysis was implemented using the Scikit-Learn library in Python and was not intended as a standalone predictive classifier, but rather as a gene-prioritization framework. The response variable was defined according to the comparison being analyzed, using either binary class labels (e.g., SA–4–1BBL vs. 3H3) or multiclass labels reflecting combinations of treatment, tissue type, and time point. Gene importance scores derived from the RF model were used to identify high-confidence regulatory candidates and were subsequently integrated with clustering results and downstream multi-omics analyses.

2.8. Integrative Analysis of Pathways, Transcription Factors, and Receptor Annotations in Immune Regulation

Pathway enrichment analysis was performed using the clusterProfiler (V4.18.4) and ReactomePA (V1.54.0) packages in R, with the enrichKEGG function used for KEGG pathways and the enrichPathway function for Reactome analysis. Adjusted p-values (Benjamini–Hochberg method) ≤ 0.05 were considered statistically significant. While gseapy was explored for accessibility and visualization in Python, the primary enrichment results presented were generated using the R-based pipeline. Transcription factors (TFs) and receptors were identified through annotation using curated databases linked to Reactome and KEGG. The enriched TFs were analyzed in the context of known biological processes relevant to the experimental design, and key TFs and receptor families were highlighted based on their annotated regulatory roles in the DEG sets rather than derived through motif or binding prediction analyses. Reactome top-level pathway proportions were calculated by categorizing all pathway annotations per cluster and determining their relative frequency (%).

2.9. IMPRes Analysis

IMPRes (Integrative MultiOmics Pathway Resolution) [44], an in-silico hypothesis generation method, was used to construct a multi-layered KEGG pathway and protein–protein interaction (PPI) network, and to overlap generated bulk RNA-seq data from SA–4–1BBL, 3H3, and Naive to uncover the most relevant sequential pathway connections for each of the studied comparisons. The IMPRes method applies a stepwise active pathway detection method using a dynamic programming approach and helps trace the connections between pathways and key genes contributing to the changes between compared datasets. The top 100-ranked DEGs for each treatment group were selected, and after removing duplicates, a total of 184 unique genes remained. To prioritize immunologically relevant genes from this refined set, the IMPRes algorithm was applied, integrating expression data with known immune-related pathways. This enrichment analysis identified 27 key genes implicated in immune activation, cytokine signaling, and co-stimulation pathways (Supplementary Figure S1). These 27 IMPRes-prioritized genes were subsequently used to perform final, treatment-specific IMPRes analyses for 3H3 and SA–4–1BBL independently. This enabled a detailed assessment of the distinct pathway activation patterns elicited by each treatment, offering insights into their unique mechanistic impacts at a systems level. The results included a pictographic tree structure illustrating the connectivity between these genes and their associated pathways, alongside a tabular output highlighting significantly enriched pathways.

3. Results and Discussion

Our analysis integrated multiple experimental parameters: treatment type, tissue source, time point, and injection number to comprehensively assess their individual and combined impact on gene expression. Rather than isolating these factors, we systematically compared 24 biologically meaningful combinations to capture shared and unique transcriptional responses, enabling a nuanced understanding of how these variables interact to shape immune modulation.

3.1. Principal Component Analysis and Differential Gene Expression Analysis of SA–4–1BBL and 3H3 Treatment

This Principal Component Analysis (PCA) represents the distribution of different samples based on their variance across PC1 (41.26%) and PC2 (18.77%), explaining a total variance of 60.03%. While PCA illustrates global transcriptional structure and overall similarity among replicates, it does not directly assess statistical power or differential expression sensitivity; therefore, differences in replicate number across experimental groups may still contribute to observed variation in DEG counts independent of biological effect size. This PCA suggests structured variability among groups, with some overlap but clear distinctions in variance contribution across principal components (Figure 1B). Differential gene expression analysis in response to 3H3 and SA–4–1BBL treatments demonstrated substantial transcriptional changes in both lymph nodes (LN) and injection sites (IS). For the 3H3 treatment, a total of 4775 DEGs were identified in LN, and 9746 DEGs in IS compared to naïve mice. After the first injection, 1970 and 36 DEGs were detected in LN on days 3 and 7, respectively, while 2429 and 715 DEGs were identified in IS at the same time points. A second injection led to substantially higher DEG counts, with 1920 and 852 in LN, and 3838 and 2764 in IS at days 3 and 7, respectively. Cumulatively, 2515 unique DEGs in response to 3H3 treatment were observed in LN and 4796 in IS, with 7311 DEGs shared between the two tissues compared to naïve (Figure 1C). These results indicate that transcriptional changes were more pronounced in IS compared to LN across both time points and injection regimens in 3H3 treatment. Moreover, the two-injection (2X) regimen was associated with both higher DEG counts and larger absolute log2 fold-change values compared to the single-injection (1X) regimen. This conclusion is based on comparative inspection of DEG tables and quantitative summaries of effect-size distributions, which revealed a shift toward higher absolute fold changes under the 2× injection condition across both tissues. In response to SA–4–1BBL, 5616 DEGs were identified in LN, and 9307 DEGs in IS compared to naïve mice. Following a single injection, 1565 and 1743 DEGs were observed at days 3 and 7 in LN, and 579 and 1045 DEGs in IS. The second injection increased DEGs to 1798 and 510 in LN and 3991 and 3692 in IS. SA–4–1BBL treatment yielded 2964 unique DEGs in LN and 5459 in IS, with a total of 8423 DEGs in the two tissues compared to naïve. The gene counts across different threshold values (Log_2_FC 0–9), showing a decreasing trend in the number of genes as the threshold increases, indicating progressive filtering based on stricter criteria.

Direct comparison analysis between SA–4–1BBL and 3H3 revealed 1453 overlapping DEGs in LN and 2057 in IS, suggesting a considerable yet tissue and time-specific transcriptional convergence. Following a single injection, 383 and 746 DEGs were detected in LN at days 3 and 7, and 515 and 303 DEGs in IS, respectively. After the second injection, DEGs in LN dropped to 108 and 216 at days 3 and 7, while IS maintained elevated levels with 486 and 753 DEGs, indicating sustained transcriptional activity. Cross-comparison analysis further identified 1153 and 1546 unique DEGs in response to SA–4–1BBL in LN and IS, respectively. As shown in Figure 1D, differential trend analysis revealed 4564 DEGs with distinct expression dynamics, including upregulated, downregulated, and oppositely regulated patterns across both treatments. Venn diagram analysis (Figure 1E) further illustrated that while substantial overlap exists between treatments, a significant number of DEGs remain uniquely regulated by either SA–4–1BBL or 3H3 treatment, and by tissue type. Notably, IS consistently showed higher DEG counts and greater divergence, underscoring its enhanced responsiveness and functional specificity to treatment. This differential trend analysis identified 4564 unique DEGs, capturing a dynamic spectrum of upregulated, downregulated, and oppositely regulated genes across tissue, treatment, and time. To clarify how transcriptional programs evolve within each treatment, we compared DEG overlap across time (D3 vs. D7) and across regimen (1X vs. 2X) separately for IS and LN. This analysis (Figure 1C) shows that differential expression across tissues and contrasts, revealing a substantially larger transcriptional perturbation at the IS than in the draining LN (total DEGs: IS = 21,110 vs. LN = 11,844). Notably, treatment-specific divergence between SA–4–1BBL and 3H3 was also more pronounced locally, with a greater number of genes showing large between-treatment differences in IS than LN (|Δlog2FC| > 1: 6502 vs. 3650), consistent with stronger agent-dependent remodeling at the IS, providing a direct view of early vs. sustained and prime vs. boost dynamics within each agent. Figure 1E summarizes DEG overlap across time (day 3 vs. day 7) and dosing (1× vs. 2× injection) for 3H3 and SA–4–1BBL. Under the 2× injection regimen, SA–4–1BBL exhibits a large shared DEG core across day 3 and day 7, whereas 3H3 shows greater condition-specific gene expression with comparatively smaller overlaps. In contrast, the 1× injection regimen is characterized by reduced DEG overlap across time for both treatments, with markedly fewer genes conserved between day 3 and day 7. Together, these patterns indicate that increased dosing is associated with enhanced persistence of transcriptional responses, with SA–4–1BBL displaying greater cross-time consistency than 3H3.

3.2. Integrated Distance Matrix and Barcode Analyses Reveal Global and Context-Dependent Molecular Signature

The distance matrix of all genes and DEGs across the 24 comparisons reveals distinct transcriptional clusters. As shown by the red-outlined clusters in the heatmap, samples within the same tissue or treatment group tend to exhibit high similarity (low-distance values), while divergent clusters, highlighted by broader matrix separation, reflect substantial transcriptional differences across time points, treatments, and injection regimens. Hierarchical clustering highlights subgroupings, particularly among immune system (IS), lymph node (LN), and naïve states, underscoring their biological relevance (Supplementary Figures S3 and S4).

To refine DEG classification by tissue, time, and treatment specifics, a 6-digit binary barcode system was developed, where each digit represents the presence or absence of differential expression in a specific condition: LN, IS, day 3, day 7, single injection, and double injection. This approach enabled the precise stratification of 4564 DEGs into 131 unique identifiers. For example, a barcode such as “1-0-1-0-1-1” would indicate that the gene is differentially expressed in LN (1), not in IS (0), at day 3 (1), not at day 7 (0), following a single injection (1), and after a double injection (1). An Upset plot (Figure 2A) was used to visualize the trend of DEGs distribution across conditions. While truncated to highlight the most prominent overlaps, it reveals that 78 barcode identifiers were shared across D3, D7, 2X, LN, and IS, encompassing 492 genes (Figure 2C, red-highlighted row in the table). Of these, 158 were more prominent in 3H3, 171 in SA–4–1BBL, and 321 were expressed in both treatments, suggesting both treatment-specific and shared transcriptional programs. Figure 2C illustrates DEG counts across specific conditions (e.g., D3_LN, D7_IS), enabling stratification by time and tissue. For instance, the “D3_2X_IS” shows 134 DEGs, with 91 upregulated in SA–4–1BBL, 28 in 3H3, and 15 showing opposite trends, demonstrating how this stratification captures biologically meaningful differences.

Figure 2D expands this into trend-based expression profiles across all 131 barcode identifiers and 4564 DEGs. These identifiers enabled precise annotation of DEG behavior across combinations of tissue, time point, and injection regimen, allowing us to detect context-dependent regulation, such as genes exclusively upregulated in IS on day 3 after a double injection. This barcode-based system offers a scalable framework for decoding how treatment, tissue, and timing collectively shape transcriptional outcomes.

3.3. K-Mean Clustering Analysis

K-means clustering of bulk RNA-seq DEGs data identified eight distinct DEGs clusters, shedding light on immune response dynamics and treatment-specific gene regulation in cancer immunoprevention. Cluster 1, with the highest gene count (1067), was predominantly influenced by the injection site (IS), while Cluster 2 (1062) reflected lymph node (LN) activity, each showing robust treatment responses. Clusters 3 and 4 revealed balanced LN and IS contributions, with Cluster 3 strongly responding to 3H3 and Cluster 4 to SA–4–1BBL. Cluster 5 displayed a unique temporal effect, evidenced by a high Day3/Day7 ratio, while Clusters 6–8, though smaller, showed specific sensitivity to time and treatment, especially with SA–4–1BBL. Color-coded analysis highlighted Clusters 1, 2, and 5 as the most responsive, with Cluster 1 leading across metrics, marking it as a central driver of immune response across conditions. This clustering analysis provides a comprehensive view of how immune responses vary with time, treatment, and location, enhancing insights into the gene regulation pathways correlated with cancer immunoprevention (Figure 3A and File S4).

Cluster 1 was particularly enriched for genes that showed a strong response to SA–4–1BBL, with key DEGs such as CD3e, CD2, CD28, Klra1, CD3e, Marco, or Tlr9, which play roles in T cell, NK cell, and innate immune system activation, supported by cytokine and transcription factor DEGs including Il18bp, Icos, or Foxp3. These genes were significantly upregulated in IS tissues, highlighting innate immune activation mediated by SA–4–1BBL. Additionally, Cd2 and Tnfsf8 genes showed tissue specificity and were more highly expressed in IS, pointing to the localized immune response generated by SA–4–1BBL Cluster 1, which was predominantly enriched in IS, featuring upregulated genes such as Tnfaip3 critical for immune regulation and cytokine production. These genes displayed elevated expression in SA–4–1BBL-treated groups, particularly in LN tissues, suggesting that LNs may play an important role in long-term immune surveillance for cancer prevention [33]. Genes involved in immune activation, such as Gata1, Crtam, and C6, are differentially expressed across conditions, emphasizing treatment-specific immune responses. Additionally, collagen-related genes (Col1a2, Loxl2) and extracellular matrix (ECM) remodeling markers indicate dynamic tissue adaptation post-treatment, particularly in LNs of SA–4–1BBL-treated mice. Metabolic shifts, including carbohydrate and lipid metabolism genes (Ckb, Galk1, Ptgds, Pparg), further distinguish the impact of these treatments. The log2 fold-change analysis indicates that SA–4–1BBL shows a more robust and prolonged activation of the immune response, particularly at later time points, such as day 7, and following two doses. This enhanced immune activation may explain why SA–4–1BBL maximum cancer immunoprevention efficacy is observed when mice are challenged with tumor cells two weeks after the last dose treatment [33] (Figure 3).

Cluster 3 includes genes enriched in both IS and LN, with notable DEGs such as Wnt9a, Asb4, and Nrep, which are associated with developmental pathways and neural immune crosstalk. Expression patterns show specificity for the 3H3 treatment, suggesting a distinct immune-stress adaptation profile. Tissue-specificity data highlights balanced immune interactions, while Mamstr and Gpr153 point to transcriptional reprogramming post-treatment. Pathway analysis shows enrichment in Wnt and signal transduction pathways, indicating a complex modulation of signaling cascades. The equal distribution across early (D3) and late (D7) timepoints and mixed 1X/2X doses suggests Cluster 3 governs transitional immune phases (Supplementary Figure S2).

Genes such as Cd8b1, Cd8a, Cd4, CD40l, Cd96, and Tcf7 dominate Cluster 4, with clear enrichment in IS tissues, pointing to enhanced cytotoxic T, CD4^+^ T cell, and NK activation and regulation, particularly after SA–4–1BBL treatment. Reactome annotations reveal significant roles in immune activation receptor signaling (Cd274) and TCR signaling, implicating these genes in activation and regulation. Treatment with SA–4–1BBL yielded the most significant upregulation, especially during the D3 timepoint and at 1X injections, indicating its role in acute effector coordination. The dynamic response underlines Cluster 4’s potential role in initiating early anti-tumor activity by enhancing and regulating immune responses.

Highly expressed genes such as H19, Upk1b, and Syt12 characterize Cluster 5, with a unique expression profile predominantly in LN tissues under SA–4–1BBL and combination treatments. These genes are linked to developmental and ECM-related functions. Pathway analysis indicates involvement in matrix remodeling and cell adhesion signaling. Notably, the expression of Gjc3 and Kcnk3 suggests alterations in cell communication and tissue permeability. Cluster 5 genes are mostly active at D3/D7 stages, supporting their role in post-activation tissue restructuring and immune priming.

Cluster 6 features genes such as Apoh, Hsd3b7, and Hibch, with marked enrichment in LN tissues and 3H3 treatment. These genes are associated with hemostasis and platelet activation pathways, suggesting an interplay between immune responses and vascular signaling. Upregulation patterns indicate 2X injections and late (D7) responses are critical. The gene expression reflects immune-vascular coordination, possibly facilitating immune cell trafficking into lymphoid tissues.

Cluster 7 is enriched in genes involved in translational control and stress responses. Rpl26, Srsf9, and Eef2 dominate this group, indicating intense cellular activity during the early immune response, particularly in IS tissues post SA–4–1BBL treatment. The cluster exhibits sharp activation at D3, hinting at involvement in initial transcriptional bursts and immune cell activation and proliferation. Reactome mapping points to ribosomal biogenesis and RNA processing. These genes highlight Cluster 7’s role in cellular readiness and activation thresholds.

Cluster 8 shows clear IS tissue dominance, with S100a8, Slpi, Ccl6, Batf3, Klrk1, Cd300a, and Lcn2 as lead DEGs, typically associated with innate immune cells and inflammation. The SA–4–1BBL treatment drives strong expression, reflecting innate immunity’s front-line role. Pathway enrichment includes innate immunity, cytokine signaling, and neutrophil degranulation, implicating the involvement of cluster in early pro-inflammatory activity. This response is strongest at D3 and under 2X conditions, indicating activation of localized defenses.

3.4. Receptor Analysis

This study demonstrates that SA–4–1BBL, a potent immune costimulatory agonist, elicits a distinct and coordinated receptor-level immune signature across tissue types, timepoints, and frequency of injection conditions. Hierarchical clustering of 312 differentially expressed receptors revealed that SA–4–1BBL induces broad transcriptional changes, particularly within IS tissue under 2X injection conditions. Cluster 1 emerged as the dominant expression module, comprising a wide array of receptors involved in NK cell and T cell activation, innate immune sensing, and inflammatory amplification in response to SA–4–1BBL treatment. This contrasted with the 3H3-treated condition, which showed narrower and more tissue-restricted profiles, most prominently represented in Cluster 3. These findings underscore the immunostimulatory capacity of SA–4–1BBL in invoking an immune response not just in the lymphoid tissues but also in the injection site (Supplementary Figure S5).

Further dissection of receptor dynamics highlighted the critical impact of both the contribution of tissues generating an immune response and the frequency of injections of the treatments on the previously published cancer immunoprevention efficacy [45]. Both tissues exhibited variable receptor expression levels dependent on the frequency and timing of treatment, suggesting differential immune engagement between peripheral and lymphoid tissues. Gene expression profiles following the 2X-injection regimen revealed a frequency-dependent pattern of receptor regulation, wherein increased treatment frequency enhanced both the breadth and magnitude of receptor activation or repression. Notably, receptors associated with innate and adaptive immune responses, such as Marco, Ccr7, Klri1, and Tnfrsf9, exhibited cluster-specific expression patterns, suggesting coordinated immune cell trafficking and inter-tissue communication over time. Collectively, this receptor-level landscape emphasizes the unique immunomodulatory function of SA–4–1BBL and supports its use as a finely tunable immunotherapeutic platform capable of driving durable and spatially coordinated immune responses (Supplementary Figure S5).

The Trend DEGs distribution analysis under Reactome pathway highlights predominant enrichment in immune system (15.6%), signal transduction (11.26%), and metabolic pathways (20.97%), emphasizing cellular communication and metabolic changes in response to treatments (Supplementary Figure S6A). The Venn diagram reveals 256 shared trend DEGs across Reactome, KEGG, and trend datasets, indicating conserved pathways, while unique distributions suggest dataset-specific functional variations. These findings reinforce the role of ECM remodeling, immune modulation, and metabolism in the analyzed system (Supplementary Figure S6B).

3.5. Enrichment of Immune and Translational Pathways

To determine the broader biological context of the transcriptional changes, we performed Reactome pathway enrichment analysis on the differentially expressed genes in Cluster 1. The dot plot (Figure 4A,B) illustrates the top enriched pathways, colored by adjusted p-value and sized by gene count. Prominent among these were interleukin-2 signaling, antigen processing-cross presentation, and DAP12 interaction, which are essential for immune activation and immune synapse formation. These pathways are consistent with increased motility and activation of myeloid cells, T and NK cells post-SA–4–1BBL treatment.

Equally notable was the enrichment of translational machinery and mRNA processing pathways, including rRNA processing, nonsense-mediated decay, and translation initiation. This supports a scenario of heightened biosynthetic demand during immune cell proliferation and effector differentiation. The upregulation of proteasome components (Psmb8, Psme1, Psmb10) and antigen-processing machinery (H2-Oa, H2-Q5, Marco) further indicates increased antigen presentation capacity associated with myeloid cell activation.

The Cluster 1 specificity heatmap (Figure 4B) reveals strong enrichment of genes involved not just innate immunity, but also in adaptive immunity, TCR signaling (e.g., Zap70, Lat, Cd3e), cytokine signaling (e.g., I2rb, Il3ra, Il18rap), interferon response genes (e.g., Ube2l6, Isg15, Stat4), and costimulation genes (e.g., Cd28, Cd40, Icos) across IS and LN samples following SA–4–1BBL treatment. These trends validate a broad activation profile impacting both innate and adaptive arms of the immune system. Also enriched were cell cycle pathways, including mitotic checkpoints and chromosome maintenance, indicating proliferative expansion of responsive immune subsets.

3.6. Pathway Enrichment Analysis

The Reactome pathway enrichment analysis and the functional pathway heatmap highlight key biological processes enriched across distinct clusters identified by K-means clustering (Supplementary Figure S2A,B). As previously noted in the clustering analysis, each cluster exhibits unique transcriptional signatures tied to specific treatments and tissue responses, emphasizing their relevance in the context of SA–4–1BBL-induced immune modulation. The Reactome enrichment results indicate a significant involvement of pathways related to extracellular matrix organization, collagen formation, immune system regulation, and muscle contraction, suggesting their biological relevance in the studied cluster (Supplementary Figure S2A). The result further supports these findings by visualizing differential expressions across key biological processes, including hemostasis, immune response, metabolism, protein modification, and muscle contraction. The red and blue color gradients represent upregulated and downregulated pathways (Supplementary Figure S2B), respectively, reinforcing the pathway enrichment trends. Notably, the presence of pathways related to platelet adhesion, immune system activation, and metabolic processes aligns with the Reactome analysis, emphasizing the interconnected roles of cellular signaling, extracellular matrix interactions, and tissue-specific functions. These findings collectively underscore the biological significance of ECM remodeling, immune modulation, and muscle contraction pathways in the dataset (Supplementary Figure S2B).

Pathway distribution analysis across clusters highlights distinct biological process enrichments, with Reactome annotations underscoring functional specializations. Cluster 1 shows significant enrichment in immune system pathways (15.65%) and extracellular matrix organization (2.34%), suggesting specialization in immune and structural functions. Cluster 2 predominantly represents metabolic pathways (26.84%) and exhibits a notable role in protein metabolism (6.40%), small-molecule transport (3.95%), and neuronal signaling (2.45%). Signal transduction and cell cycle pathways are concentrated in Cluster 6 (15.18% and 16.52%, respectively), with additional contributions to gene expression (6.70%) and vesicle-mediated transport (4.46%). Developmental biology is associated with Cluster 8 (2.82%). Cluster 5, though less broad in overall enrichment, shows distinct engagement with the muscle contraction pathway (15.5%) and DNA repair (4.5%), suggesting a possible link to tissue-specific stress responses or local immune–tissue interactions. Comparatively, Cluster 2 is metabolically oriented (26.8%), with significant activity in protein metabolism (6.4%) and small molecule transport (4.0%), while Cluster 6 is enriched for signal transduction (15.2%), cell cycle regulation (16.5%), and gene expression (6.7%), indicative of proliferative and transcriptionally dynamic states.

These findings echo the earlier K-means clustering results, where each cluster demonstrated distinct transcriptional patterns shaped by treatment modality and tissue specificity. The Reactome-based functional profiling thus reinforces the biological significance of the identified clusters, revealing complex inter-cluster relationships and underscoring the immune-metabolic diversity orchestrated by SA–4–1BBL treatment (Figure 4C).

3.7. In Silico Hypothesis Generation Using IMPRes

While traditional enrichment analysis identifies immune pathway involvement, it does not reveal how these pathways interact or which key genes orchestrate the transitions between them. To address this limitation and uncover mechanistic cascades, we employed IMPRes using the differentially expressed genes (DEGs) from 3H3 and SA–4–1BBL conditions as input. IMPRes was run separately for each condition to reconstruct pathway-level networks and identify key connectors bridging innate and adaptive immune pathways (Figure 5A,B). The IMPRes output enabled directional mapping of immune signaling hierarchies, providing context for how specific genes such as Stat1, Cd247, Ifng (in SA–4–1BBL), and Casp1, Ccr5, Cxcl10 (in 3H3) function as pivotal mediators of their respective immune programs.

The IMPRes results from gene expression-based comparisons between 3H3 and SA–4–1BBL treatments reveal distinct immune-modulatory effects relevant to cancer immunoprevention. SA–4–1BBL-treated samples exhibited significantly higher expression of key chemokines (Ccl7, Ccl4) and cytokines (Il10, Tnf, Il1b), suggesting enhanced immune cell recruitment and activation (Supplementary Figure S8). Temporal analysis showed sustained pro-inflammatory responses in SA–4–1BBL conditions, with elevated Cxcl10 and Tnf levels persisting through day 7. In contrast, 3H3 conditions showed comparatively dampened expression of inflammatory and metabolic regulators, such as Mapk12 and Ppara. These results highlight the potential of SA–4–1BBL to enhance anti-tumor immunity by activating chemotactic and pro-inflammatory pathways (Supplementary Figure S8).

Pathway enrichment analysis revealed that SA–4–1BBL activated pathways involved in Th1 and Th2 cell differentiation, T cell receptor signaling, and NK cell-mediated cytotoxicity. Genes such as Ifng, Cd247, and Stat1 showed significantly higher expression in SA–4–1BBL treatment, consistent with the IMPRes-inferred network highlighting these as key connectors downstream of TCR activation and cytokine signaling. Additionally, genes including the regulator of interaction between T cells and dendritic cells (Tnfrsf11a, Il3ra), and inflammatory bowel disease or T cell differentiation markers (Il4, Il12rb1) were enriched, indicating broader immunomodulatory roles for SA–4–1BBL. These effects were further supported by heatmaps showing elevated expression of adaptive immune genes (Figure 5D).

In contrast, 3H3 treatment induced expression patterns associated with innate immunity, including Toll-like receptor signaling (Stat1, Cxcl10), NOD-like receptor signaling (Casp1, Il1b), and chemokine signaling (Ccr5, Cxcr2). IMPRes revealed that these genes formed a tightly connected inflammatory cascade centered on Calm4, Camk2d, and Casp1, suggesting a skew toward inflammation and potential neurodegenerative signatures (Figure 5A). Notably, genes implicated in tuberculosis and neurodegeneration (e.g., Camk2d, Calm4) were more highly expressed in the 3H3 group, indicating a distinct immunological orientation.

Heatmaps of log-fold changes across conditions (Figure 5C,D) further confirmed distinct gene expression patterns, where adaptive immune response genes (Ifng, Cd247, Il4, Tbx21, Stat1) were upregulated under SA–4–1BBL treatment, whereas innate immune and pro-inflammatory genes (Casp1, Ccr5, Il1b) were enriched in the 3H3 condition. Interestingly, genes involved in NK cell cytotoxicity (Cd244, Sh2d1b1, Vav1) were upregulated in both conditions but with consistently higher levels in SA–4–1BBL-treated samples, suggesting a more potent cytotoxic immune profile.

Hierarchical clustering confirmed these trends, showing that SA–4–1BBL-treated samples formed a distinct transcriptional cluster associated with adaptive immune activation, while 3H3-treated samples clustered with signatures of innate immune stimulation and inflammation. These patterns, taken together with IMPRes-reconstructed signaling cascades, provide mechanistic insight into how SA–4–1BBL drives T cell immunity, whereas 3H3 shapes a more inflammatory microenvironment (Figure 5A–E).

A Venn diagram comparison of the two IMPRes outputs (Figure 5E) revealed limited pathway overlap, with only three shared pathways (chemokine signaling, cytokine-cytokine receptor interaction, NK cell-mediated cytotoxicity), and a larger set of uniquely enriched pathways for SA–4–1BBL (e.g., T cell receptor signaling, Th17 cell differentiation) and 3H3 (e.g., NOD-like receptor signaling, cAMP signaling). These results underscore the distinct immune modules engaged by each treatment.

Overall, these analyses reveal a compartmentalized immune architecture, where gene clusters and pathway modules respond uniquely based on both tissue context (IS vs. LN) and treatment modality. SA–4–1BBL drives a prolonged, transcriptionally driven immune program enriched in tissue remodeling, adaptive immunity, and innate immunity, as suggested in enrichment analysis, while 3H3 activates proinflammatory cascades associated with innate immunity.

3.8. Cancer Immunoprevention Resource (CIR) Web Portal

To facilitate further exploration and validation of our findings, we have deposited the analyzed bulk RNA-seq data in the KBCommons database [45], an open-access resource for multiomics datasets. Researchers can access the analyzed data, including differentially expressed genes, pathway enrichment results, and heat maps, through the KBCommons portal (https://kbcommons.org:8443/?organism=MusMusculus, accessed on 8 January 2025). This repository allows users to query gene expression levels, compare pathway activation across conditions, and integrate these results with other publicly available datasets. By making our data accessible, we aim to support further investigations into the cancer immunoprevention mechanisms of SA–4–1BBL and encourage future studies, including single-cell analyses, to build upon our findings. The Cancer Immuno-Prevention Resource web portal (https://cir.missouri.edu) was created to visualize the plots based on experimental design and bioinformatics analysis (Figure 6).

4. Conclusions

This study was conducted to provide insight into the mechanisms underlying the immunoprevention efficacy of SA–4–1BBL through integrated transcriptomic profiling. We have previously reported that treatment of naïve mice with 100 μg of SA–4–1BBL as a single agent, twice two weeks apart, followed by tumor challenge two weeks later, showed the most robust cancer immunoprevention efficacy, ranging from 60% to 100% efficacy without detectable toxicity. In contrast, treatment with an equimolar dose of 3H3 agonistic Ab failed to protect mice from tumor development and, in a side-by-side study with SA–4–1BBL, showed significant toxicity. Limited mechanistic studies demonstrated the critical role of IFN-γ and CD4^+^ T and NK cells in the observed cancer immunoprevention. To further elucidate the cellular and molecular players that contribute to SA–4–1BBL-mediated immunoprevention and to develop this molecule for clinical translation, we used the same cancer immunoprevention treatment for RNA-seq analysis of the injection site tissue and draining LNs to generate bioinformatics data for a more comprehensive insight into immune changes and formulate a testable hypothesis. A central strength lies in the spatially resolved transcriptomic analysis of LNs and IS tissues to uncover their complementary roles in orchestrating immune responses.

K-means clustering identified eight gene expression clusters, with Clusters 1 and 4 dominant in IS tissues and enriched in effector-related genes such as Cd4, Cd40, Cd40l, Cd28, Klra1, Lat, and Cd8b1. These genes drive cytotoxic cell, CD4^+^ T, and NK cell activation, immune synapse formation, and pro-inflammatory signaling. Elevated levels of Ccl4, Ccl7, and Tnf further support IS as a hub of early immune activation and effector recruitment to promote an innate-to-adaptive cascade. Conversely, LN-specific clusters (e.g., Cluster 2) were enriched in immune stimulatory, regulatory, and tissue-remodeling genes, including Tnfaip3, Irf5, Crtam, and Col1a2, that potentially generate a more balanced immune surveillance response without collateral damage. The expression of genes such as H19 and Syt12 suggests enhanced immune retention and long-term memory development within the LN microenvironment, which is consistent with a prolonged 12–14 weeks of immune protection window against cancer in our published studies.

The integrated analysis of cluster-wise boxplots and statistical summaries reveals clear, condition-dependent transcriptional responses. Cluster 8 stands out with a high mean Log2 Fold Change (~0.52) and significant p-values across all conditions, reinforcing its role as a consistently upregulated cluster under both control and immune-stimulated conditions, possibly representing genes involved in core inflammatory or activation pathways. Conversely, Cluster 2 and Cluster 5, with mean Log2 Fold Changes around −1.4 and −3.2, respectively, show significant downregulation, particularly in early immune stimulus (D3_3H3_LN) and across SA–4–1BBL comparisons, suggesting strong transcriptional repression or immune evasion-related activity. Cluster 7, with large interquartile ranges and condition-specific spikes (e.g., strong upregulation in D7_SA–4–1BBL_IS), reflects highly stimulus-specific gene activation, likely tied to late immune or costimulatory responses (Supplementary Figure S7). Importantly, statistical significance is universal across clusters (all ANOVA p-values = 0), validating the observed visual trends. Significant q-values for comparisons such as D7_SA–4–1BBL_vs_D7_3H3_IS across Clusters 1–8 suggest robust differential gene expression in response to SA–4–1BBL stimuli. Cluster 3 and Cluster 4 also demonstrate variable responses, but their median Log2FC values (~0.14 and ~0.006) and moderate significance suggest nuanced modulation rather than overt activation or repression. In summary, this combined evidence supports a model where certain clusters (e.g., 1, 8, 6, 7) act as positive responders to immune signaling, while others (e.g., 2, 5, 3) are negatively regulated or repressed, forming a bidirectional transcriptomic shift. These findings lay a foundation for functional annotation and potential biomarker or therapeutic target discovery.

Temporal profiling showed that IS responses peaked early (Day 3), aligning with initial stimulus encounter, whereas LN responses were more prominent at later stages (day 7 and post-boost), consistent with memory responses and regulation. Reactome pathway analysis reinforced these dynamics: IS clusters were enriched in TCR signaling, innate cytokine networks, and T and NK cell activation pathways, while LN clusters exhibited signatures of adaptive memory immune regulation, including protein synthesis, interleukin signaling, and metabolic programming. Differential gene expression between SA–4–1BBL and 3H3 revealed 1453 DEGs in LN and 2057 in IS, with immune-stimulatory genes such as Tnfrsf14, Cd28, Cd40l, Cd40, Il12b, and Il2ra more highly expressed in SA–4–1BBL-treated tissues. This indicates the communication between innate and adaptive immune responses. Key mediators such as CD40l and Klra1 enhance stimulation and T cell activation, while Tnfaip3 functions as a regulatory hub, preventing overactivation and systemic inflammation (Supplementary Figure S8). Importantly, SA–4–1BBL selectively induces chemokines (Ccl2, Ccl3, Ccl5) and cytokines (Tnf, Il10, Il12b, Il18rap), promoting immune recruitment without triggering systemic toxicity. Transcriptional comparison confirmed that 3H3 drives rapid activation of innate pathways, C-type lectin receptors, Toll-like receptor signaling, and neurodegeneration-associated networks, suggesting a role in acute inflammatory response. In contrast, SA–4–1BBL promotes adaptive regulation through Th1–Th2 differentiation, cancer-related pathways, apoptosis, and inflammatory disease signaling. Despite some shared signatures, such as cytokine–receptor interaction and NK cell-mediated cytotoxicity, SA–4–1BBL produced a more sustained immunological effect, evidenced by elevated Cxcl10 and enriched T cell receptor signaling. Heatmap analysis further demonstrated discrete gene expression patterns, underscoring SA–4–1BBL’s role in promoting myeloid cell activation, Th1/Th2 differentiation, and tumor-targeted immunity. These observations are consistent with our previous studies, which demonstrated significant dysregulation of immune responses and a systemic cytokine storm generated by 3H3 Ab, but not in SA–4–1BBL.

Our integrative transcriptomic profiling confirms that gene expression clusters align with treatment type, tissue localization, and immune function. Differences in pathway enrichment and interaction types revealed distinct immunological trajectories, such as transient receptor modulation (3H3) versus sustained signaling, modulatory activation, and proliferation (SA–4–1BBL). Clusters enriched in Cd28, Tnfaip3, Crtam, and Irf5 correlated with expanded T cell and NK cell populations, while increased expression of ECM genes (Col1a2, Loxl2) paralleled observed tissue remodeling and fibroblast activation. IS-associated Clusters 1 and 4 (SA–4–1BBL) mapped to cytotoxic cells, CD4^+^ T and NK cells, while LN-dominant Cluster 2 (SA–4–1BBL) overlapped with metabolically active myeloid cells. This mapping confirms the cellular specificity of bulk-defined clusters and validates their interpretability. A limitation of this study is the lack of independent validation of protein expression levels for selected DEGs, and the reliance on PCA-based clustering, which reflects global expression structure but does not fully account for replicate-dependent differences in differential gene detection sensitivity; these aspects will be addressed in future studies through targeted protein-level validation and expanded experimental designs.

Collectively, the transcriptional programs identified here align with and expand on the previously reported limited mechanisms underlying the immunoprevention efficacy of SA-4-1BL. The transcriptomic analyses define a compartmentalized, spatiotemporal immune architecture induced by SA–4–1BBL, where the injection site appears to serve as a hub for immune activation and coordinates with the LN to sustain and balance the immunosurveillance response. The controlled, durable, and tissue-specific immune stimulation underscores the potential of SA–4–1BBL as a safe and effective cancer immunoprevention agent and highlights the utility of system-level approaches for dissecting the complexity of immunopreventive responses. Overall, these bioinformatic analyses provide a hypothesis-generating framework that highlights candidate pathways underlying immunoprevention efficacy of SA–4–1BBL. This will structure future studies aimed at mechanistic and functional validation.