Abstract
Insufficient infiltration of cytotoxic T cells into solid tumors
remains a critical obstacle in cancer immunotherapy. Despite extensive
efforts to comprehend the mechanisms governing this limited
infiltration, few studies have focused on the evolution of T cell
motility behavior after co-culture. In this study, we combined
quantitative cell trajectory analysis, computational modeling, and
bulk/single-cell RNA sequencing to systematically characterize the
impact of cell interactions. We reveal that in a 2.5D co-culture system
with multiple cancer-cell clusters, cancer-specific T cells exhibit
increased directional persistence, which facilitates their efficient
searching of cancer-cell clusters. Additionally, these T cells form
prolonged interactions with cancer cells, which is the most crucial
factor for their accumulation on cancer-cell clusters. Furthermore,
post-interaction, a cancer-cell subpopulation displays
immunosuppressive traits, reducing T cell attractant expression, and
undergoing epithelial-to-mesenchymal transition. These findings offer
valuable insights into improving immunotherapy efficacy and tackling T
cell infiltration challenges in solid tumors.
Subject terms: Systems analysis, Computational models, Computer
modelling, Tumour immunology
__________________________________________________________________
Quantitative analysis of T-cell trajectories and gene expression in
cancer and T cells reveals how cancer-T cell interactions shape the
spatiotemporal dynamics of antigen-specific T-cell motility and drive
tumor immune evasion in co-culture.
Introduction
Cytotoxic T lymphocytes (CTLs, CD8^+) have emerged as a pivotal force
in cancer immunotherapy due to their remarkable cytotoxic
capabilities^[38]1,[39]2. However, challenges persist in optimizing the
effectiveness of immunotherapies^[40]3–[41]5. For instance, cancer
cells can exploit mutations to evade immune surveillance, thus
diminishing the efficacy of immunotherapy^[42]6,[43]7. Previous studies
have underscored the impact of extrinsic components, such as
tumor-associated macrophages^[44]8,[45]9, myeloid-derived suppressor
cells^[46]10,[47]11, cancer-associated fibroblasts^[48]12,[49]13, and
the extracellular matrix^[50]14,[51]15, in hindering the infiltration
of CTLs in solid tumors, thereby reducing patient response to
immunotherapy^[52]16,[53]17. Despite these findings, few studies have
delved into intrinsic factors, such as the motility behaviors and gene
expression profiles of CTLs with different tumor-recognition abilities.
These observations highlight the necessity of investigating CTL
motility within the tumor microenvironment. While recent studies have
emphasized the essential role of collective behaviors of CTLs in tumor
killing^[54]18,[55]19, several questions remain unanswered. For
instance, what strategies do CTLs employ when navigating multi-tumor
clusters, balancing between exploration and exploitation^[56]20?
Additionally, what are the temporal changes in gene expression among
CTLs and cancer cells and how do they further affect the behavior and
interactions of these cells? Answers to these questions will be crucial
for enhancing T-cell engineering and design to address obstacles in
cancer immunotherapy^[57]21,[58]22.
Technically, it is challenging to characterize cellular behaviors in
vivo, particularly due to the necessity of intravital high-resolution
real-time imaging techniques and the constraints on imaging
duration^[59]23,[60]24. To address these challenges, several 3D in
vitro co-culture systems, such as spheroids, scaffold-based models, or
microfluidics platforms have been developed to investigate the
interactions between T cells and tumors^[61]18,[62]25,[63]26, as well
as T cells and extracellular matrix fibers^[64]14,[65]27,[66]28. These
systems serve as valuable tools for studying T cell mobility patterns
and their interactions within controlled environments. While 3D systems
offer a more physiologically relevant environment, they often require
more sophisticated imaging equipment and considerable technical
expertise^[67]29,[68]30. Consequently, 2.5D systems have been developed
to study collective cellular behavior^[69]30,[70]31, despite certain
apparent limitations in fully replicating the physical and biochemical
microenvironment of real tissues. Nevertheless, 2.5D systems could
provide a viable alternative, requiring more accessible imaging
technology and facilitating detailed investigations of cellular
behaviors.
In this study, we harnessed a 2.5D in vitro co-culture platform to
quantify the influence of cancer-T cell interactions on the
spatiotemporal dynamics of CTL motility. Using an integrated approach
combining live-cell microscopy, quantitative image analysis,
computational modeling, and transcriptomic profiling, we decoded the
temporal motility patterns of CTLs in proximity to and atop multiple
cancer-cell clusters post interaction. Furthermore, we quantified the
impact of these spatiotemporal dynamics on CTL’s cancer-seeking
behavior and their accumulation on cancer-cell clusters. In addition,
gene expression analysis revealed potential molecular drivers
underlying CTL behavioral adaptations and diverse tumor-immune
responses. Our research underscores the critical role of cancer-T cell
interactions in governing CTL movement patterns, offering valuable
insights into the spatial organization of CTLs encircling cancer cell
clusters within solid tumors.
Results
A time-lapse pipeline for both motility and gene-expression of CTLs
To investigate the time-dependent interactions between cancer cells and
cancer-specific CTLs and their impact on CTL mobility patterns, spatial
distribution, as well as gene-expression, we developed a comprehensive
pipeline, integrating time-lapse imaging, cell tracking and trajectory
analysis, computational modeling based on cell mobility statistics,
time-resolved gene expression profiling of both cancer cells and CTLs
(Fig. [71]1A).
Fig. 1. Overview of the co-culture system and pipelines.
[72]Fig. 1
[73]Open in a new tab
A Schematic representation illustrating the integrated workflow for
experimental set-up, motility analysis, and gene expression profiling.
B Diagrams of the co-culture experiments and representative 2D
microscopic images for the WT and OVA groups. Scale bars, 100 µm. C
Quantification of the number of recruited CTLs on cancer cells over
time for WT or OVA groups, manually counted within a 15 µm surrounding
region every 30 min (see Methods). For each time point, the total
number of CTLs residing on cancer cells within the imaging field was
counted (n = 28 for the OVA groups and n = 9 for the WT groups from one
independent experiment, mean ± s.e.m.). D Counts of live cancer cells
over time in different regions of the imaging field for the WT and OVA
groups (n = 9, mean ± s.e.m.). E Interaction time of CTLs on cancer
clusters for B16-OVA and B16-WT groups (n = 217 vs. 225). C–E
Statistical significance was determined using the Two-sided Wilcoxon
rank-sum test: p < 0.001 (***), p < 0.01 (**), p < 0.05 (*), p > 0.05
(n.s.). For (C), the p values at t = 0 h, 6 h, 12 h, and 18 h were
2.13E-01, 8.26E-02, 4.96E-05 and 3.72E-04, respectively. For (D), the p
values at t = 0 h, 6 h, 12 h, and 18 h were 3.32E-02, 5.22E-03,
4.11E-05 and 4.11E-05, respectively. For (E), the p-value is 4.46E-41.
We employed a co-culture model in a 2.5D µ-slide chamber, combining
pre-activated OT-1 CD8^+ T cells with either OVA-expressing B16
melanoma cells (B16-F10-OVA) or wild-type B16 cells (B16-F10-WT) as
controls (Fig. [74]1A, B, Methods). After embedding the CTLs in a thin
collagen matrix overlaying a sub-confluent tumor monolayer
(Supplementary Fig. [75]1A), time-lapse imaging was performed for 24 h
at 2-min intervals, allowing continuous observation of cell dynamics
within a ~5.27 µm optical section. Our interest lies in understanding
the behavior of T cells within such thin layer of collagen-fiber
network around cancer-cell clusters (Supplementary Fig. [76]1A), and
the following time-lapse imaging analysis is specifically tailored to
track their movement around and on tumor cells within the focal imaging
plane.
Qualitatively, under the microscope, the most apparent observation was
cancer cells could be recognized by CTLs labeled with the red
fluorescence in Fig. [77]1B (Supplementary Fig. [78]1B, Supplementary
Movie [79]1 and Supplementary Movie [80]2). Not surprisingly, the
density of CTLs increased more rapidly over time in co-culture with OVA
cancer cells than in co-culture with WT cancer cells (Fig. [81]1C and
Supplementary Data [82]1). In addition, the WT cancer cells
demonstrated continued growth in the co-culture with CTLs, whereas the
number of OVA cancer cells gradually decreased (Fig. [83]1D and
Supplementary Data [84]1).
To elucidate the faster accumulation of CTLs shown in Fig. [85]1C, we
first hypothesized that it might result from the faster growth or
swifter mobility of CTLs co-culture with OVA cancer cells. However, no
significant differences in CTL doubling time or migration velocity were
observed in the two co-culture conditions (Supplementary Fig. [86]1C,
Supplementary Fig. [87]1D and Supplementary Data [88]1). Alternatively,
the accumulation of CTLs could be attributed to their reduced mobility
around the OVA cancer cells or a more pronounced chemotactic motion
towards OVA cancer cells. Intriguingly, manually tracking CTLs on both
types of cancer cells revealed an extended engagement between CTLs and
OVA cancer cells (Fig. [89]1E). Given the constraints of manually
handling CTL tracks, in the following analysis, we turned to an
automated tracking strategy to investigate the potential chemotaxis
behavior of CTLs.
Quantitative motility patterns of antigen-specific T cells
To quantitatively investigate the motility behaviors of CTLs in the
co-culture with OVA and WT cancer cells, we focused on characterizing
their motility patterns before and after CTLs reached cancer cell
clusters. The workflow for time-lapse image processing and data
analysis is illustrated in Fig. [90]2A, B (Methods, Supplementary
Fig. [91]2). A selection criterion was applied to include only T cells
exhibiting continuous motion across at least 10 frames, thus
guaranteeing the consistency and validity of motion analysis within the
imaging plane. Utilizing the automated image analysis pipeline, we
obtained a total of 4276 and 5146 CTL tracks from five biological
replicates in the OVA and WT groups, respectively. Subsequent analyses
were conducted based on these CTL trajectories and the spatial
boundaries of cancer-cell clusters (Methods).
Fig. 2. Quantitative motility patterns of CTLs in the co-culture with OVA and
WT cancer cells.
[92]Fig. 2
[93]Open in a new tab
A Pipeline for processing time-lapse images and tracking T cells. B
Schematic representation of data analysis, including dwell time, mean
square displacement, directional bias, directional persistence of CTLs
towards cancer cell clusters. C Density of CTLs on cancer clusters
(cells/mm^2) over time, demonstrating a significant T cell accumulation
in OVA groups compared to WT groups (n = 5). Statistical significance
was determined using the Two-sided Wilcoxon rank-sum test: p < 0.001
(***), p < 0.01 (**), p < 0.05 (*), p > 0.05 (n.s.). The p-values at
0 h, 6 h, 12 h, and 18 h were 1, 1, 7.94E-03 and 7.94E-03,
respectively. D Dwell time distribution of T cells on cancer-cell
clusters in Series 1, divided into 6-h intervals. The left panel
(0–6 h) shows n = 149 adhesion events for the WT groups and n = 110
events for the OVA groups; the right panel (6–12 h) shows n = 131
events for the WT groups and n = 223 events for the OVA groups. The
distribution for 12–18 h is shown in Supplementary Fig. 3D. E Step size
distributions of CTLs over 6-hour intervals: Left panel (0–6 h), middle
panel (6–12 h), and right panel (12–18 h). The p-values for (E) (left,
middle, and right) are 9.9E-23, 3.3E-26 and 1.8E-70, respectively. Data
represent results from n = 5 independent experiments. D, E Statistical
significance was determined using the Two-sided two-sample
Kolmogorov–Smirnov test. KS statistic D and p-value are shown inside
the plotting boxes. F Directional bias and G persistence distributions
of CTLs movement towards the nearest cancer cell boundary within three
time-intervals. Statistical significance was determined using the
two-sided unpaired Student’s t test with p < 0.001 (***), p < 0.01
(**), p < 0.05 (*). [a, b) defines the value range where x ≥ a and
x < b. C–G Red lines (or bars) represent OVA groups, and cyan lines (or
bars) represent WT groups. C, F, G, n = 5 independent experiments,
mean±s.e.m. The p-values for (F, G) are provided in Supplementary
Data [94]8.
Firstly, we quantified the behaviors of CTLs when their fluorescent
signals overlapped with those of cancer cells, indicating direct
physical contact and interactions between CTLs and cancer cells
(Methods). A marked increase in CTL accumulation on OVA cancer-cell
clusters was observed within 6–12 h post co-culture (Fig. [95]2C,
Supplementary Fig. [96]3A, B and Supplementary Data [97]2), consistent
with earlier manual counting (Fig. [98]1C). Notably, the total number
of CTLs per unit area remained similar between the two co-culture
systems over time (Supplementary Fig. [99]3B), suggesting that enhanced
accumulation of CTLs on OVA cancer-cell clusters is unlikely due to
differential CTL proliferation, in agreement with prior observations.
Although we occasionally observed CTLs “hopping” onto cancer clusters
from other z-planes, quantification revealed that such events were
rare—approximately 30 times less frequent than those approaching cancer
cells within the thin fiber network layer (Supplementary Fig. [100]3C).
Furthermore, to evaluate CTL-cancer cell interaction dynamics, we
quantified the dwell time of CTLs on cancer-cell clusters
(Fig. [101]2D; Methods). The distribution closely followed a
heavy-tailed power-law pattern. Notably, within the first 6 h, dwell
time distributions of CTLs were similar between co-cultures with OVA
and WT cancer cells (Kolmogorov-Smirnov test, D = 0.15, p = 0.096).
However, during 6–18 h of the co-culture, CTLs in the OVA condition
exhibited a noticeably longer-tailed distribution (D = 0.25,
p = 6.2 × 10^−5), with some CTLs remaining on cancer-cell clusters for
up to 100 minutes (Fig. [102]2D, Supplementary Fig. [103]3D–H and
Supplementary Data [104]2). These findings suggest that once CTLs
recognized and interacted with cancer cells, they tended to engage for
prolonged periods.
Secondly, to investigate the motility behaviors of CTLs before they
reached cancer-cell clusters, we studied the step sizes, mean square
displacement (MSD), directional bias and directional persistence of
CTLs whose fluorescent signals did not overlap with cancer cells
(Methods). We observed that the difference in step-size distribution
between the two groups gradually increased over time, with the KS
statistic rising from D = 0.05 at 0–6 h to D = 0.07 during 12–18 h, as
determined by the two-sample Kolmogorov–Smirnov test (p < 0.001,
Fig. [105]2E). After 12 h, CTLs co-cultured with OVA cancer clusters
exhibited shorter step sizes compared to those with WT cancer clusters,
particularly within 5 μm of the cancer boundary (D = 0.001, p < 0.001,
Supplementary Fig. [106]4). Moreover, CTLs demonstrated super-diffusive
motion, characterized by
[MATH: MSD~τα :MATH]
, with the exponent α around 1.3 (Supplementary Fig. [107]5A).
Intriguingly, during the first 6 h of co-culture, no significant
difference in α was observed between the two co-cultures (Supplementary
Fig. [108]5A, panel 2). However, in the subsequent 6–12 h (panel 3) and
12–18 h (panel 4), CTLs co-cultured with OVA cancer cells exhibited
higher α values compared to those with WT cancer cells, suggesting that
CTLs co-cultured with OVA cancer cells may possess a more efficient
searching capability^[109]32.
To further investigate the motility patterns of CTLs with respect to
their distance towards cancer cell clusters, we examined directional
bias and directional persistence (Methods, ref. ^[110]33). CTLs
co-cultured with OVA cancer cells gradually increased their motion bias
towards the nearest cancer-cell clusters from the first 6 h to the
9–18 h of the co-culture, particularly within 0–10 μm of the
cancer-cell boundary (Fig. [111]2F, Supplementary Fig. [112]5B and
Supplementary Data [113]2). By contrast, CTLs co-cultured with WT
cancer cells showed a gradual decrease in directional bias, especially
within 0–10 μm range. Although CTLs co-cultured with OVA cancer cells
only displayed a mild bias (~0.1) toward the nearby cancer clusters
within 9–18 h (Fig. [114]2F and Supplementary Fig. [115]5B), the
difference was significant compared to CTLs co-cultured with WT cancer
cells.
For the directional persistence, CTLs co-cultured with OVA cancer cells
exhibited considerably stronger persistence compared to those cultured
with WT cancer cells, particularly during 3–12 h (Fig. [116]2G,
Supplementary Fig. [117]5C and Supplementary Data [118]2). These
observations indicate that the interactions between OVA cancer cells
and CTLs might enhance the ability of CTLs to maintain the consistent
movement direction. It is noteworthy that a consistent decrease in
directional persistence over time was observed in both co-culture
systems. One possible explanation is that the orientation or density of
collagen fibers surrounding cancer cell clusters might affect the
directional persistence. However, we did not observe significant
differences in fiber orientation around cancer cell clusters between
the two co-culture conditions (Supplementary Fig. [119]6A–C).
Therefore, the mechanism underlying such decrease of persistent motion
when CTLs approaching cancer-cell clusters is still elusive.
Furthermore, we assessed pairwise correlations among step size,
directional persistence, and directional bias across varying distances
from the cancer boundary per 6-hour intervals (Supplementary
Fig. [120]7). No significant correlation was found between step size
and directional bias. Although mild correlations (~0.2) between
directional persistence and bias were observed within 10 μm of cancer
clusters, they were negligible. Moreover, step size and directional
persistence exhibited similarly weak correlations (~0.25) beyond 10 μm
in both co-culture conditions. Overall, these findings suggest that the
interplay among these motility parameters is weak and
distance-dependent.
Computational studies revealed the key motility patterns governing efficient
searching and the accumulation of CTLs on cancer-cell clusters
In Fig. [121]2, we provided a quantitative and systematic
characterization of the 2-dimentional motility patterns of CTLs. To
further elucidate the roles of different T-cell motility patterns in
their searching and accumulation on cancer-cell clusters, we developed
a 2-dimensional agent-based model where four motility characteristics
of T cells, i.e., dwell time, directional bias, persistence, and
step-size distribution, were varied and their effects on the
accumulation of T cells were compared (Fig. [122]3A; Methods), with a
display of simulated trajectories in one simulation in Fig. [123]3B and
Supplementary Movie [124]3.
Fig. 3. Computational studies revealed the key motility patterns for the
efficient searching and the accumulation of CTLs on cancer-cell clusters.
[125]Fig. 3
[126]Open in a new tab
A Schematic representation of agent-based models. B An illustration
depicting the trajectories of 20 CTLs after moving 100 steps each. C
Schematics of the calculation method of mean arrival time of CTLs
towards cancer clusters. D Violin plot and statistical significance
test of the mean arrival time of CTLs to cancer clusters under
different patterns when altering only the differences in the
directional persistence or directional bias or step sizes. 100
replicated simulations were performed, with each simulation consisting
of 200 cells moving for 800 steps. Statistical significance was
determined using the two-sided unpaired Student’s t-test: p < 0.001
(***), p < 0.01 (**), p < 0.05 (*), p > 0.05 (n.s.). The p-values from
left to right are 2.97E-05, 1.59E-05 and 4.36E-02, respectively. The
accumulation of CTLs on cancer-cell clusters as a function of time when
considering different variations: E step size, F directional
persistence, G directional bias, and H dwell time. The shaded area
represents the 95% confidence interval of the simulated data, while the
bold line represents the mean. Statistical significance was determined
using the Two-sided Wilcoxon rank-sum test: p < 0.001 (***), p < 0.01
(**), p < 0.05 (*), p > 0.05 (n.s.). For clarity in visualization, data
points are displayed with sub-sampling at an interval of 5. For (E),
the p-values at t = 100, 200, 400, and 600 were 1.77E-3, 3.31E-2, 0.651
and 0.398, respectively. For (F), the p-values at t = 100, 200, 400,
and 600 were 0.315, 0.113, 0.651, and 0.398, respectively. For (G), the
p-values at t = 100, 200, 400, and 600 were 2.97E-06, 3.26E-12,
1.16E-08 and 1.04E-06, respectively. For (H), the p-values at t = 100,
200, 400, and 600 were 2.27E-26, 4.98E-25, 1.33E-22 and 2.20E-24,
respectively.
Firstly, we explored the impact of motility patterns on the mean
arrival time of CTLs towards cancer-cell clusters (Fig. [127]3C).
Specifically, the distance-dependent directional persistence, bias, and
step-size distribution of CTLs characterized from the two coculture
systems were utilized (Methods). By manipulating one pattern while
maintaining the other two patterns consistent with the WT co-culture
system, we investigated whether the observed differences in CTL
motility patterns between the two co-culture systems could influence
the mean arrival time of CTLs. Notably, statistically significant
differences were observed when altering the directional persistence or
bias, rather than the step-size (Fig. [128]3D and Supplementary
Data [129]3). Intriguingly, CTLs in the co-culture with OVA cancer
cells reached cancer-cell clusters more rapidly, emphasizing the
importance of directional persistence and bias in CTLs’ search
strategy.
In addition, to ensure the robustness of our findings, we conducted
similar analyses with a single cancer-cell cluster in the simulation
(Supplementary Fig. [130]8A, B & Supplementary Movie [131]4). The
observed trends remained consistent when considering only the
differences in directional persistence or step-size, while showing
inconsistency in directional bias, which can be attributed to the
initial distribution of CTLs (Supplementary Fig. [132]8C). These
results suggest that stronger directional persistence could enhance
CTLs’ efficiency in searching for cancer-cell clusters.
Subsequently, we investigated the impact of motility patterns on the
accumulation of CTLs on cancer-cell clusters. No significant difference
in the number of CTLs on cancer-cell clusters was found when varying
the step-size (Fig. [133]3E) or the directional persistence
distribution (Fig. [134]3F), with dwell-time distribution held constant
as in the WT co-culture. In contrast, altering directional bias
resulted in a modest change in CTL accumulation (Fig. [135]3G and
Supplementary Data [136]3). Notably, simulations incorporating the
dwell-time distribution from the OVA co-culture resulted in a marked
increase in CTL accumulation, compared to those using the WT-derived
distribution. (Fig. [137]3H). This effect remained qualitatively robust
regardless of cancer-cell clusters configuration (Supplementary
Fig. [138]8D–G). These findings suggest that enhanced adhesion, rather
than the motility outside of cancer-cell clusters, plays a dominant
role in the accumulation of CTLs on cancer clusters.
Bulk RNA sequencing analysis revealed potential molecular mechanisms
underlying T-cell motility patterns
To investigate the potential molecular mechanisms underlying T cell
motility patterns, we performed both bulk and single-cell RNA
sequencing of cancer and T cells sorted from the co-culture systems
(Methods). The PCA plot following batch effect correction is presented
in Supplementary Fig. [139]9A. Each point represents a biological
replicate collected from multiple microwells across at least two
independent experiments (see Methods for details). Note that the two
C0A-R1 samples (B16-F10-OVA cells before co-culture, Round 1; hereafter
referred to as “R1”) are further away from the C0-R1 (B16-F10-WT cells
before co-culture, R1), C0-R2 (B16-F10-WT cells, Round 2; “R2”), and
C0A-R2 (B16-F10-OVA cells, R2) samples. This separation may indicate an
inconsistency in the sequencing results of the C0A-R1 samples, which
are further addressed in the Discussion part.
First, analysis of bulk RNA-seq data of T cells co-cultured with
B16-F10-OVA cancer cells revealed that, in addition to classical
activation and cytotoxicity-associated cytokines such as interleukin-2,
perforin, interferons, and TNFs, chemokines potentially involved in the
recruitment and polarization of innate immune cells—such as Xcl1, Ccl1,
Ccl3, Ccl4, and Csf2—were also significantly upregulated
(Fig. [140]4A), showing a gradual induction after co-culture initiation
(Fig. [141]4B). Overall, 262 genes were upregulated at 12 h post
co-culture (Supplementary Data [142]4). KEGG and Gene Ontology (GO)
enrichment analyses of these genes revealed significant enrichment in
cytokine-cytokine receptor interaction and chemokine signaling pathways
(Fig. [143]4C, D), suggesting that antigen-specific CTLs may
orchestrate the local immune microenvironment through recruitment of
innate immune components.
Fig. 4. Bulk RNA sequencing analysis revealed potential molecular mechanisms
underlying T-cell motility patterns.
[144]Fig. 4
[145]Open in a new tab
A Volcano plot of differentially expressed genes between CTLs
co-cultured with OVA cancer cells and those co-cultured with WT cancer
cells at 12 h after the co-culture (T12A vs T12). Genes that were
significantly up and down regulated in T12A samples, i.e., log2(fold
change) > 1 and FDR < 0.05, were shown in red and blue, respectively. B
Temporal dynamics of selected up-regulated genes in (A). CTLs at 0, 6,
12 h of the co-culture with OVA or WT cancer cells were filtered and
sequenced in two replicates. The Z-score of the corresponding
expression of the selected genes were calculated and demonstrated. C
Results of the KEGG pathway enrichment analysis of differentially
expressed genes shown in (A), where the size of the circle corresponds
to the number genes in the specific pathway and the color corresponds
to FDR values of the enrichment analysis. D Results of the gene
ontology (GO) enrichment analysis of differentially expressed genes
shown in (A), where different colors correspond to GO terms for
biological process (BP), molecular function (MF) or cellular component
(CC). E Volcano plot of differentially expressed genes between OVA and
WT cancer cells at 12 h after the co-culture with CTLs (C12A vs C12).
Genes that were significantly up and down regulated in C12A samples
were shown in red and blue, respectively. F Temporal dynamics of
selected up-regulated genes in (E). OVA and WT cancer cells at 0, 6,
and 12 h of the co-culture with CTLs were filtered and sequenced in two
replicates. The Z-score of the corresponding expression of the selected
genes were calculated and demonstrated. G Results of the KEGG pathway
enrichment analysis of differentially expressed genes shown in (E). H
Results of the gene ontology (GO) enrichment analysis of differentially
expressed genes shown in (E).
For the gene expression profiles of cancer cells, compared to
B16-F10-WT ones, 942 genes were upregulated in B16-F10-OVA cancer cells
at 12 h after the co-culture with CTLs (Fig. [146]4E, Supplementary
Data [147]4). Given the upregulation of interferon-gamma by
antigen-specific CTLs, it is not surprising that many of these
genes—including Stat3, Ifitm3, and Cdkn1a—are associated with
IFN-responsive pathways. To assess whether these differences were
intrinsic to the two tumor cell lines or induced by CTL interaction, we
performed additional analyses. First, we exposed B16-F10-OVA and
B16-F10-WT cells to IFN-γ (50 ng/mL, 24 h) and analyzed their
transcriptomes. The two cell lines exhibited highly overlapping gene
expression changes (Supplementary Fig. [148]9B), with no significant
differences in pathways related to cell cycle, antigen presentation,
immune regulation, or apoptosis (Supplementary Fig. [149]9C). In
addition, EdU incorporation assays showed comparable proliferative
activity between the two lines under baseline conditions (Supplementary
Fig. [150]9D). Furthermore, flow cytometry analysis of MHC-I surface
expression demonstrated similar levels between the two cell lines at
baseline (Supplementary Fig. [151]9E). Together, these results suggest
that the transcriptional divergence observed in co-culture divergence
observed in co-culture could be largely driven by antigen-specific CTL
interactions, rather than inherent differences in IFN-γ responsiveness
or baseline functional states between the tumor cell lines.
Notably, the chemokine genes Cxcl9 and Cxcl10, critical for T cell
recruitment, were progressively upregulated in OVA cancer cells
following CTL interaction, with expression levels markedly higher than
those in the CTLs themselves (Fig. [152]4F). Moreover, in response to
CTL-mediated attack, cancer cells upregulated genes enriched in
pathways such as antigen presentation, apoptosis, and TNF signaling,
along with GO terms like MHC complex and cell killing (Fig. [153]4G, H,
Supplementary Fig. [154]9F). Consistently, apoptosis-related genes were
significantly elevated in B16-OVA cells after co-culture, compared to
baseline (Supplementary Fig. [155]10), reflecting a strong CTL-induced
cytotoxic response.
Furthermore, we noted that pathways of cell adhesion molecules and
ECM-receptor interaction are enriched for genes up-regulated in OVA
cancer cells (Supplementary Fig. [156]9G). Notably, several genes
related to epithelial-to-mesenchymal transition (EMT) were upregulated,
such as Serpine1, Col1a1, Dab2, etc. Using GSEA with human homologs
from the 50 MSigDB hallmark gene sets, we observed a strong enrichment
of the EMT signature (EMT, NES:1.77, FDR:5.6e-05, Supplementary
Fig. [157]9H). Consistently, EMT-related GO terms were also
significantly enriched (Supplementary Data [158]4). These observations
suggest that in response to T-cell attack, cancer cells may tend to
metastasize. To evaluate differences in migratory behavior between
B16-F10 and B16-F10-OVA cells, we performed a wound healing assay
before and after co-culture with CTLs (Supplementary Fig. [159]11). In
the absence of co-culture, both cell lines exhibited comparable
migration. Notably, following co-culture, B16-F10-OVA cells displayed
significantly enhanced wound closure, indicating an increase in
migration capacity. In contrast, B16-F10 cells showed reduced
migration. These observations suggest that co-culture with CTLs
differentially affects the migratory behavior of the two cell lines,
although further studies using dedicated metastasis assays would be
required to assess their metastatic potential (Supplementary
Fig. [160]11).
Single-cell RNA sequencing analysis unveiled heterogeneity and resistance
mechanisms
By analyzing the bulk RNA-seq data of co-cultured cancer and CTLs, we
have revealed essential genes that were differentially expressed in OVA
cancer cells and antigen-specific CTLs, respectively. On the other
hand, we also observed from the time-lapse images that some cancer
cells survived longer and some CTLs remained engaged with cancer cells
for extended periods. To further investigate the molecular basis of
this heterogeneity within each cell type, we performed single-cell RNA
sequencing (scRNA-seq) on the sorted cancer cells and CTLs.
First, we applied the scRNA-seq clustering of all cells at a resolution
of 0.5, and then projected the results on a UMAP (Methods), identifying
nine sub-clusters of cells (C1–C9, Fig. [161]5A). Prior to sequencing,
CD8^+ T cells were separated from tumor cells using Dynabeads (11462D,
ThermoFisher), allowing us to categorize the cells into four groups:
g-B16-F10, g-CTLs-B16-F10 (CTLs co-cultured with B16-F10),
g-B16-F10-OVA, and g-CTLs-B16-F10-OVA (CTLs co-cultured with
B16-F10-OVA). These groups are visualized from left to right in the
UMAP panels (Fig. [162]5A). Based on gene signatures, cells in clusters
C1 and C2 were categorized as B16-F10 cancer cells after 12 h of
co-culture with CTLs, while cells in clusters C3 and C4 were identified
as B16-F10-OVA cancer cells after 12 h of co-culture with CTLs.
Clusters C5–C9 were assigned as CTLs. It is important to note that, due
to the ability of T cells to attach to cancer cells, sorted cell groups
may include overlapping clusters. Specifically, CTLs in clusters C5,
C6, and C7 were identified within the sorted group g-B16-F10, while
CTLs in clusters C5, C6, C7, and C9 were identified within the sorted
group g-B16-F10-OVA. Interestingly, a distinct cluster of CTLs (C9) was
observed in association with sorted OVA cancer cells but not with WT
ones (Fig. [163]5A).
Fig. 5. Single-cell RNA sequencing analysis unveiled heterogeneity and
resistance mechanisms.
[164]Fig. 5
[165]Open in a new tab
A UMAP of scRNA-seq data of the sorted cancer and T cells. All cells
were combined to perform the clustering, and 9 cell clusters (C1 to C9)
were identified after removing ones with low RNA contents or high
mitochondrial RNAs (Methods). WTcancer and WTctl (as well as OVAcancer
and OVActl) cells were isolated after 12 h of co-culture using
Dynabeads positive selection with EpCam and CD8 antibodies,
respectively. B Top ten differential expressed genes of the cell
clusters shown in (A). The color represents the average expression of
cells in the indicated cell cluster. The size of the circle represents
the percentage of cells that express the DEG in the indicated cell
cluster (Methods). C Merge of non-unique CTL-subclusters and rename of
cancer and CTL clusters. OVActl cells belonging to CTL-subclusters that
overlapped with the ones in WTctls were merged into one new subcluster
named as CD8T-cx. D Gene-expression differences between OVActl cells
that belong to CD8T-c2 and CD8T-cx in (C). Genes expressed in more
cells (x-axis) from CD8T-c2 and CD8T-cx as well as a higher average
expression (y-axis) were shown in red and blue, respectively. E
Gene-expression differences between Cancer-c3 and Cancer-c4 cells in
(C). Genes expressed in more cells (x-axis) from Cancer-c3 and
Cancer-c4 as well as a higher average expression (y-axis) were shown in
red and blue, respectively. F Ligand-receptor cellular interactions
between CTLs (CD8T-c2) and cancer cells (Cancer-c4 or -c3) in the
coculture of OVA cancer cells and CTLs (Methods). Only the interactions
with a p value smaller than 0.01 were displayed. The color represents
the communication probability as indicated by the color bar.
Following the identification of the nine clusters, we conducted
differential expression gene (DEGs) analysis and generated a heatmap to
visualize transcriptional differences across clusters (Methods and
Supplementary Fig. [166]12A). For cancer cells, WT clusters (C1–C2) and
OVA clusters (C3–C4) exhibited distinct expression patterns, consistent
with bulk RNA-seq results. For CTLs, clusters C5–C8 showed highly
similar profiles, while cluster C9 displayed notable differences,
indicating a distinct subpopulation within OVA-CTLs. Consistently, the
top 10 signature genes of the subcluster C9 (Fig. [167]5B) were among
the up-regulated genes identified in bulk RNA-seq (Fig. [168]4A), but
exhibited significantly lower expression in the other OVA-CTL
sub-clusters. This result suggested that only a subset of OVA-CTLs
(mainly C9) contributed significantly to the observed differences from
WT-CTLs in bulk analysis. In addition, volcano plot analysis between
CTLs in C9 (CD8T-c2 in Fig. [169]5C) and the merged cluster of C5–C8
(CD8T-cx) revealed that CD8T-c2 cells exhibited elevated expression of
chemokine genes (e.g., Xcl1, Ccl1, Ccl3, and Ccl4), adhesion-related
genes (e.g., Crtam) and effector genes (e.g., Ifng and Prf1), with a
higher proportion of cells expressing these genes (Fig. [170]5D and
Supplementary Data [171]5). Functional enrichment analysis of CD8T-c2
DEGs highlighted pathways involved in positive regulation of cytokine
production, cell killing, and regulation of cell–cell adhesion
(Supplementary Fig. [172]12B), in agreement with bulk RNA-seq analysis.
For cancer cells, scRNA-seq analysis also revealed heterogeneity within
OVA-cancer cells, particularly in response to CTLs attack. Notably, the
Cancer-c4 subcluster exhibited elevated expression of migration- and
metastasis-associated genes, such as Fn1^[173]34, Tubb3^[174]35, and
Dst^[175]36, compared to Cancer-c3 (Supplementary Data [176]5). Among
them, Fn1 was expressed in a 15% higher proportion of Cancer-c4 cells.
Although Tubb3 and Dst showed similar percentages of expressing cells
between the two clusters, their average expression levels were
approximately two-fold higher in Cancer-c4 (Fig. [177]5E). Conversely,
Cancer-c3 cells showed higher expression of T-cell chemokine genes
(Cxcl9, Cxcl10) and antigen-presentation-related genes (H2-T23, H2-K1)
than those in Cancer-c4. For instance, Cxcl9 exhibited approximately
two-fold higher expression in Cancer-c3, despite a similar proportion
of cells expressing the gene.
Functional enrichment analysis further supported these transcriptional
differences. DEGs in Cancer-c3 were enriched in pathways associated
with interferon response, positive regulation of cytokine production,
and cytokine-mediated signaling response (Supplementary Fig. [178]12C).
By contrast, DEGs in Cancer-c4 were enriched in pathways were enriched
in pathways related to cell adhesion, ameboid cell migration,
regulation of cell morphogenesis, etc. (Supplementary Fig. [179]12D,
Supplementary Data [180]5). Gene set enrichment further confirmed a
strong enrichment for EMT-related genes in Cancer-c4 (Supplementary
Fig. [181]12E). These results, together with the wound healing assay
(Supplementary Fig. [182]11) showing enhanced migration in B16-F10-OVA
(Cancer-c3/c4) compared to B16-F10 (Cancer-c1/c2), suggest that
Cancer-c4 cells may acquire higher metastatic potential and CTL
resistance upon interaction.
Finally, to move beyond transcriptional changes and explore the
intercellular communication underlying cancer-T cell interactions, we
performed CellChat analysis^[183]37 analysis using the scRNA-seq data
of cancer and T cells at 12 h after the co-culture. Notably, Cancer-c3
cells have several unique interactions with the CTLs (Fig. [184]5F),
such as CD6-ALCAM, H2-K1-CD8A/CD8B1 and H2-T23-CD8A/CD8B1. Some of
these interactions have been shown to be important for the
activation/expansion/adhesion of T cells^[185]38,[186]39 or the
antigen-presentation of cancer-cells^[187]40,[188]41. Moreover, it was
interesting to note that only CTLs in CD8T-cx cluster interacted with
OVA cancer cells via CXCL9-CXCR3 and CXCL10-CXCR3, which could be
responsible for the motility behavior of antigen-specific T cells
toward cancer cells (Supplementary Fig. [189]12F).
To validate the role of the CXCR3-ligand axis (comprising CXCL9,
CXCL10, and CXCL11) in regulating CTL motility, we conducted co-culture
experiments of B16-F10-OVA tumor cells and CTLs using the CXCR3
antagonist ACT-660602 (HY-151096, MCE)^[190]42. Specifically, 408 nM
ACT-660602 was added to the co-culture medium, while the control group
was maintained under identical conditions without the antagonist
(wi/wo-ACT-OVA). Following ACT-660602 treatment, CTLs in the co-culture
system exhibited reduced motility persistence compared to untreated
controls (Supplementary Fig. [191]13A, B) suggesting a critical role
for the CXCR3-ligand axis in guiding CTL migration toward tumor
targets. Notably, no significant differences were observed in dwell
time between the treated and control groups (Supplementary
Fig. [192]13C), indicating that CXCR3 signaling primarily modulates
locomotory behavior rather than cell-cell interaction dynamics.
Furthermore, a quantitative analysis of the temporal dynamics of tumor
cell area and T cell density on tumor-cell clusters revealed that the
ACT-660602-treated group exhibited a larger tumor cell area and reduced
T cell density on tumor-cell clusters at 18 h post-co-culture compared
to the untreated group (Supplementary Fig. [193]13D–F). These findings
suggest that the inhibition of the CXCR3-ligand axis impairs both CTL
recruitment efficiency and subsequent effector-mediated tumor
elimination.
Discussion
In this work, we combined a 2.5D co-culture system with an automated
quantitative analysis pipeline to characterize the motility patterns of
cytotoxic T cells co-cultured with cancer cells expressing a specific
antigen. We found that antigen-specific CTLs exhibited greater
directional persistence outside of cancer-cell clusters and prolonged
dwell time on cancer cells compared to non-specific CTLs. In addition,
through computational simulations, we showed that increased directional
persistence can expedite the first arrival time of CTLs on cancer
cells, while longer dwell time can augment the accumulation of CTLs on
cancer-cell clusters. These motility patterns of CTLs might facilitate
their localization and coordinated assault on cancer cells.
To decipher the potential mechanisms underlying the observed motility
patterns of CTLs, we also performed gene expression analysis and
founded that both cancer cells and CTLs expressed more chemokines and
cytokines after the attack of CTLs. These molecules might increase the
directional persistence and bias of CTLs towards cancer-cell clusters.
Further knock-out experiments are warranted to quantify the effects of
these molecules on CTL motility. Previous studies have shown that
CCL3/4 is responsible for attracting CTLs in 3D co-culture
systems^[194]19. In this work, by blocking the CXCR3-CXCL9/10
signaling, the results suggested that CXCL9/10 could also play a
crucial role in the observed enhanced directional bias of CTLs in our
2.5D co-culture system. Note that, in addition to the chemotaxis-driven
biased motion, at a distance close to tumor-cell clusters, a higher
directional bias of CTLs could also result from CTLs’ direct contact
interactions with tumor cells (Supplementary Fig. [195]14), potentially
triggering TCR-mediated changes in motility^[196]19,[197]43. Therefore,
it would be valuable to investigate this phenomenon further using
systematic gene knockout experiments and higher-resolution imaging
techniques^[198]44,[199]45 in the future.
Furthermore, to understand the extended dwell-time of CTLs on cancer
cells, we identified several potential adhesion candidates mediating
interactions between cancer cells and antigen-specific CTLs, such as
the MHC-I/TCR complex and integrins, based on gene expression data.
Notably, TCR activation has been reported to reduce the T cell motility
via various mechanisms^[200]46–[201]48. The quantitative contribution
of these molecular pathways on dwell time could be further consolidated
using the gene knock-out, RNAi or interaction-blockade experiments in
future studies.
In addition to CTLs, our analysis of scRNA-seq data revealed that
cancer cells differentiated into two groups following interaction with
antigen-specific CTLs. One group exhibited reduced expression of
antigen-presenting molecules and chemokines but increased expression of
EMT-related and immunosuppressive genes. These observations might help
to provide mechanistic insights on the potential immune evasion of
cancer cells. Previous studies have suggested that cancer-cell EMT
induced by immune attack can suppress the activity of
CTLs^[202]49,[203]50. In parallel, recent studies on sublethal T cell
interactions highlight the possibility that tumor cells can acquire
resistance following non-lethal contacts with T cells^[204]51,[205]52.
Together, detailed mechanistic studies are still required to unveil the
detailed biological mechanisms. Additionally, a limitation is that not
all biological replicates for bulk RNA-seq were processed in fully
independent runs, and the single-cell RNA seq was performed only once.
Nevertheless, the consistency of our findings is supported by IFN-γ
response comparisons (Supplementary Fig. [206]9B), DEG overlap analysis
(Supplementary Fig. [207]9F), and cross-validation between bulk and
single-cell datasets (Supplementary Data [208]7). This limitation could
be further addressed in future studies involving 3D co-culture systems.
Finally, our investigation focused solely on the interactions between
cancer cells and antigen-specific CTLs in the controlled 2.5D
co-culture. However, actual tumors comprise various other cell types,
such as cancer-associated fibroblasts, endothelial cells, and
tumor-associated macrophages. Co-culturing multiple cell types in vitro
holds promise for elucidating the effects of paracrine and autocrine
interactions on the spatial distribution of antigen-specific
CTLs^[209]53. In addition, the current computational model remains a
simplification of the complex in vivo environment. For instance, it
does not take into account cell proliferation, cytotoxic killing, or
tumor evolution processes such as metastasis. Moreover, the influence
of fiber network distributions—potentially affecting both motility and
cell-cell interactions—was also not incorporated. Future models that
integrate these factors may offer a more comprehensive view of CTL
dynamics within tumor microenvironments. While certain motility
behaviors may differ in vivo (Supplementary Figs. [210]15–[211]16), the
pathway-level similarity in molecular responses suggests that the
candidate molecules mediating T cell navigation and sustained adhesion
present promising targets for T cell engineering, potentially enhancing
the success of cancer immunotherapy. Additionally, the procedure to
quantify motility patterns of T cells in the 2.5D co-culture system
might be useful for studying T-cell motility behavior in 3D co-culture
systems.
Methods
Cell samples
B16-F10 murine melanoma cells were kindly provided by the Stem Cell
Bank, Chinese Academy of Sciences (see Supplementary Data [212]6). To
label the cancer cells, B16-F10 cells were first transduced with a
plasmid encoding the Ca^2+ sensor GCaMP6s
(pLV[Exp]-Puro-EF1A > GCaMP6s) using lentiviral infection.
Subsequently, B16-OVA cells were generated by transducing these cells
with a plasmid encoding ovalbumin (pLV[Exp]-Neo-EF1A > {ovalbumin}) via
lentiviral infection. The transduced B16 cancer cells were then
cultured in DMEM (GIBCO) supplemented with 10% FBS (GIBCO) and 1%
penicillin/streptomycin (GIBCO).
OT-1 CD8^+ T cells isolation, and activation
CD8⁺ T cells were isolated from the spleens of female C57BL/6J-OT-I
mice (Cyagen Biosciences, Suzhou; Cat# C001198). Single-cell
suspensions were prepared using Cell Staining Buffer, and CD8⁺ OT-1 T
cells were positively selected using CD8 magnetic beads. The isolated
cells were then activated using the Dynabeads™ T-Activator CD3/CD28 kit
(GIBCO, Cat#11456D) according to the manufacturer’s protocol. After
activation, the beads were removed by magnetic separation and washing
to obtain purified activated T cells.
The naive CD8⁺ OT-I T cells were activated using Dynabeads® Mouse
T-Activator CD3/CD28 (GIBCO)^[213]54 according to the manufacturer’s
instructions. Cells were cultured in RPMI 1640 medium supplemented with
heat inactivated 10% FBS (GIBCO), 1% penicillin/streptomycin (GIBCO),
2 mM L-glutamine (GIBCO), 50 mM β-mercaptoethanol (GIBCO) and 30 U/ml
mouse recombinant IL-2 (Biolegend).
Flow cytometry assay of OT-1 CD8^+ T cells
For flow cytometry analysis, the activated OT-1 CD8^+ T cells were
fixed and stained with fluorophore-conjugated antibodies diluted in
flow cytometry staining buffer. After staining, the cells were
resuspended in buffer and immediately acquired using a CytoFLEX S flow
cytometer (Beckman). Data were analyzed using FlowJo software. (1)
Lymphocytes were initially gated based on forward scatter (FSC) and
side scatter (SSC) parameters. (2) Single cells were selected within
the lymphocyte gate to exclude doublets and cell aggregates. (3) CD3
expression was used to identify the T cell population. (4) The CD3⁺
population was further gated to identify CD8⁺ T cell subsets.
Quantitative analysis indicated that CD8⁺ T cells accounted for ~98.8%
of the population, confirming the high purity of the isolated T cells
(Supplementary Fig. [214]17).
Co-culture of cancer and T cells
Co-culture experiments were conducted in 2.5D µ-slide chambers (ibidi,
8-well), using pre-activated, OVA-specific CD8⁺ cytotoxic T lymphocytes
(OT-1 CTLs) and B16-F10 mouse melanoma cells. Two tumor cell lines were
used: B16-OVA (expressing Ovalbumin) and B16-WT (lacking OVA
expression) as the control. B16-F10-OVA or B16-F10-WT cells were seeded
onto Poly-D-lysine (Sigma) coated µ-Slide wells and incubated
overnight. For antigen presentation, B16-OVA cells were pulsed with
1 μM OVA peptide at 37 °C for 1 h, then washed twice with serum-free
DMEM. Next, OT-1 CTLs were labeled with anti-CD45-PE (Invitrogen) and
mixed with neutralized collagen I solution (Corning) at a final E:T
(effector-to-target) ratio of 2:1. The collagen solution was prepared
by adjusting the pH to 7.5 using NaOH and mixing with 10× DPBS. The
CTL-collagen suspension was added onto the sub-confluent tumor
monolayer. After polymerization for 45 min at 37 °C, 400 µL of complete
medium was added to each well. The experimental setup was adapted from
protocols described in ref. ^[215]55.
Live-cell imaging-based cytotoxicity assay
Time-lapse imaging was performed for 24 h at 2-min intervals. Tumor
cells were visualized via the GCaMP6s fluorescent reporter, while CTLs
were labeled with anti-CD45-PE through immunofluorescence staining.
Bright-field and time-lapse fluorescent images were acquired using a
NIKON TI2-E system equipped with a 20× objective lens. FITC
(excitation: 480 nm; emission: 512–550 nm) and TxRed (excitation:
540 nm; emission: 577–632 nm) channels were used for fluorescence
detection. All co-culture experiments were conducted at 37 °C and 5%
CO[2] within an incubation chamber enclosing both the microscope stage
and body.
Sample preparation for bulk and single-cell RNA sequencing
After co-culture for 0 h, 6 h, and 12 h, a collagenase solution was
added to the µ-Slide 8 well and incubated for 15 min to release T cells
from the collagen gel, followed by trypsin EDTA for 5 min to detach
cancer cells. The resulting cell suspension containing cancer cells and
T cells was subjected to positive selection using Dynabeads conjugated
with EpCam and CD8 antibodies, respectively. Bead-bound cells were then
resuspended in release buffer and the bead-free cell suspension was
transferred to new tubes. Bulk RNA sequencing was performed on cells
from co-culture at 0 h, 6 h, and 12 h, whereas single-cell RNA
sequencing was specifically conducted on cells obtained from the 12 h
co-culture.
For bulk RNA-seq, samples were derived from multiple biological
replicates across at least two independent experiments. In each
experiment, approximately 12,000 tumor cells were seeded per well in
µ-Slide 8-well plates for cytotoxicity assays. To obtain sufficient
material (~500,000 cells per condition per time point), we cultured
over 10 plates in parallel for each condition (OVA and WT at 0, 6, and
12 h), referred to as the first batch (“R1”). This procedure was
independently repeated to generate a second batch of samples (“R2”), to
assess the transcriptional responses of B16-F10 and B16-OVA cell lines
to IFN-γ stimulation (C0-R2, C0A-R2, WT-IFN, OVA-IFN). For single-cell
RNA-seq, samples were collected from one experiment, using pooled cells
from ~40 wells per condition.
EdU incorporation assay
To evaluate potential differences in proliferative activity between the
two cancer cell lines, a cell cycle analysis was performed using an EdU
incorporation assay. B16-F10 and B16-OVA cells were seeded onto 35-mm
glass-bottom dishes and incubated for 24 h. EdU labeling was conducted
using the VF555 Click-iT EdU Universal Cell Proliferation Detection Kit
(MCE, Cat#HY-K1086), following the manufacturer’s instructions.
Fluorescent signals were imaged using a Nikon Ti2-E fluorescence
microscope.
MHC class I surface expression Assay
To assess MHC class I surface expression levels in B16-F10 and B16-OVA
cell lines, flow cytometry analysis was performed using a PE-conjugated
anti-mouse H-2Kd antibody (Biolegend, Cat#116607). Cells were stained
according to the manufacturer’s instructions and analyzed on a
FongCyte™3 Flow Cytometer. Data acquisition and subsequent analysis
were carried out using FlowJo software.
Wound healing assay
B16-F10-OVA cells post-CTL killing (clusters c3 and c4), B16-F10-WT
cells post-CTL killing (clusters c1 and c2), and untreated B16-F10-WT
and B16-F10-OVA cells (as controls) were seeded in 96-well plates at a
density of 3000 cells per well and cultured for 24 h. Once the cells
reached confluency, a straight scratch was introduced into the
monolayer using a sterile cell scraper. The wells were then washed with
PBS to remove debris and replenished with fresh medium. Time-lapse
bright-field images of the wound were captured every 10 min over a 24-h
period using a NIKON Ti2-E microscope system equipped with a 20 ×
objective. The wound area was quantified by manually tracing the
cell-free regions in sequential images using ImageJ software, and
closure rates were calculated as the percentage reduction in wound area
over time.
Validation of CXCR3-ligand axis functionality
To investigate CXCR3-ligand axis functionality, we implemented targeted
inhibition using ACT-660602 (HY-151096; MCE), a selective CXCR3
antagonist. Specifically, 408 nM ACT-660602 was supplemented into the
culture medium of the B16-F10-OVA co-culture system (wi-ACT-OVA), while
the control group maintained identical conditions without the
antagonist (wo-ACT-OVA). Immune-mediated cellular interactions were
continuously tracked through time-lapse imaging following standardized
protocols from our live-cell imaging-based cytotoxicity assay.
Mouse experiments
All experiments were conducted in accordance with the ethical
guidelines approved by the Animal Research Committee of Guangdong
Medical University (Protocol No. GDY2104034). All mice were maintained
under SPF (Specific Pathogen-Free) conditions. Six-week-old female
NOD-SCID mice were purchased from Cyagen Biosciences (Suzhou) Inc
(Cat#C001070). To evaluate the anti-tumor properties of OT-1 CD8 + T
cells activated with anti-CD3/anti-CD28 Dynabeads in vivo, melanoma
mouse models were established using NOD-SCID mice. A total of 20 mice
were injected with either B16-F10-OVA (GCaMP6s-transfected) (n = 10) or
B16-F10 (GCaMP6s-transfected) cells (n = 10), with the latter serving
as the control group. Cage assignments were not randomized to minimize
animal stress from frequent handling; however, all treatments and
measurements were performed within consistent daily time windows to
reduce potential confounders.
Tumor model establishment
For melanoma model induction, approximately 5 × 10^5 B16-F10-OVA
(GCaMP6s-transfected) (n = 5) or B16-F10 (GCaMP6s-transfected) cells
(n = 5) were suspended in 50 µL PBS, mixed with an equal volume of
Matrigel (50 µL), and injected subcutaneously into the dorsal aspect of
the right hind leg of mice. When tumor volumes reached approximately
100 mm³ (approximately 7 days post-injection), each tumor-bearing mouse
received 5 × 10^6 OT-1 CD8 + T cells via tail vein injection. Body
weight, tumor dimensions and general health status were monitored every
four days. Humane endpoints included tumor volume >1 cm³, weight loss
>20%, ulceration, or distress signs. Mice meeting these criteria were
euthanized immediately; others were euthanized 15 days post-injection
using approved methods. To ensure objective quantification, tumor
volume measurements were performed independently by two trained
technicians blinded to group allocation. All raw measurement data were
compiled and anonymized by a third researcher uninvolved in
experimental procedures. One B16-F10 melanoma-bearing tumors exceeded
the ethical maximum tolerated volume (1 cm³) and was excluded from
analysis. Statistical analysis was performed using the two-sided
unpaired Students t-test (n ≥ 4).
Cell isolation for comparative transcription analysis
To compare T cell gene signatures between in vitro co-culture systems
and in vivo tumors, B16-F10-OVA (n = 5) and B16-F10 (n = 5)
melanoma-bearing mice were euthanized two days post-T cell
administration. Fresh tumor tissues were minced into ~1 mm fragments
using sterile scalpels and enzymatically dissociated in digestion
medium containing mouse collagenase (1 mg/mL) and hyaluronidase
(1000 U/mL) at 37 °C for 45 min. The resulting cell suspension was
filtered through Falcon 100 µm cell strainers to remove undigested
tissue, washed repeatedly to eliminate dead cells, and resuspended in
culture medium (DMEM supplemented with 10% FBS, 1%
penicillin-streptomycin, 1% sodium pyruvate, 1% HEPES, 500 µL insulin,
and 10 nM estrogen). Cell concentration and viability were quantified
using Countess™ Cell Counting Chamber Slides.
OT-1 CD8 + T cells were positively selected from the suspension using
CD8+ antibody-conjugated Dynabeads. After magnetic separation,
bead-bound cells were incubated with release buffer, and bead-free
CD8 + T cells were collected for bulk RNA sequencing. Concurrently,
GCaMP6s-expressing B16-F10-OVA and B16-F10 tumor cells were isolated
via fluorescence-activated cell sorting (FACS) and processed for bulk
RNA sequencing.
Manual analysis of CTL recruitment, division time and interaction time
For measurement of effector cell recruitment (Fig. [216]1C), the
detected regions of cancer cells were expanded by adding a belt area
with a constant diameter of 15 μm (“dilation circle”), and CTLs within
the dilated region were manually counted. Division time and interaction
time of OT-1 CTLs were analyzed by Fiji^[217]56,[218]57 (Supplementary
Fig. [219]1C and Fig. [220]1E). Migration velocity was processed using
Imaris software (RRID:SCR_007370) (Supplementary Fig. [221]1D).
Quantification of CTL trajectories
CTLs were identified using the TxRed fluorescence channel and their
trajectories were generated utilizing Trackmate 7^[222]58, an
open-source Fiji plugin. First, the T cell spots were detected using
the LOG (Laplacian of Gaussian) detector. The estimated object diameter
was set to 14 μm and the quality threshold was set to 0.18.
Subsequently, track segments were established using the LAP tracker,
with the maximal distance allowed for frame-to-frame linking set to
25 μm. Additionally, only tracks with a duration of at least 10 frames
were selected for further analysis. Finally, the resulting track table
and spot table were obtained. The former provides information about the
trajectories of the CTLs over time while the latter contains the
coordinates of all T-cell centroids.
Identification of cancer-cell clusters
Cancer-cell clusters were imaged by splitting the FITC fluorescence
channel. To extract the main tumor region, basic image processing was
conducted using Fiji. Tumor cell positions were updated frame-by-frame
to ensure that their dynamic changes were accurately incorporated into
the analysis. This encompassed threshold adjustment, Gaussian blur, and
background subtraction. Subsequently, the preprocessed tumor images
underwent further processing using the MATLAB Image Processing Toolbox
and customized algorithms. These operations included morphological
closing, hole filling, edge smoothing, and the removal of small
impurities. Cancer-cell clusters were defined as connected regions with
an area exceeding 400 pixels. These operations were then applied to
stacks of tumor images, generating a time sequence of binarized
cancer-cell clusters data.
Quantification of the density of CTLs
To quantify the accumulation of CTLs on cancer-cell clusters over time
(Fig. [223]2C), the density of CTLs on cancer-cell clusters was
calculated by dividing the number of CTLs present on cancer-cell
clusters by the area of tumor regions, measured in cells/mm^2. The
positions of CTLs obtained from Trackmate tracking relative to
cancer-cell clusters were determined by assessing whether the centroids
of CTLs overlapped with tumor regions.
The density of CTLs around the cancer clusters (cells/mm^2,
Supplementary Fig. [224]3A) was calculated by dividing the number of T
cells within each dilation belt (enlarged to a diameter of 5 μm) by the
area of the dilation belt, which expanded every 5 μm from the boundary
of the cancer clusters, covering a total range of 0–45 μm.
Additionally, the density of T cells on cancer clusters was also
included. A 3D representation of CTL density around cancer clusters was
displayed for 6 h, 9 h and 12 h, respectively.
The total density of CTLs within the imaging field (Supplementary
Fig. [225]3B) was quantified by dividing the total number of CTLs in
the imaging field by the total area of the field, also measured in
cells/mm^2.
Dwell-time distribution of CTLs on cancer-cell clusters
To quantify the duration of CTL interactions with OVA cancer cell and
WT cancer cell, the dwell-time distribution of CTLs was assessed
(Fig. [226]2D, Supplementary Fig. [227]3D–H). Firstly, a dwelling event
was defined as the occurrence when the centroid of CTL was within a
minimal distance of less than 5 μm from the boundary of cancer-cell
clusters. This criterion was established based on radius of the CTL is
approximately 5–6 μm. The dwell-time was then calculated by determining
the number of consecutive frames in which CTL was present on
cancer-cell clusters. In cases where a CTL attached to the tumor cell
clusters multiple times within a track, each attachment was considered
as a separate dependent adhesion event.
To compare the dwell time distributions between CTLs cocultured with
OVA cancer cell and those with WT cancer cell within different stages,
the adhesive time distributions during three specific time periods were
statistically calculated: 0–6 h, 6–12 h, and 12–18 h.
Directional bias and persistence analysis
To explore the motility patterns of CTLs with respect to their distance
towards cancer cell clusters in the two co-culture systems, we
conducted a study utilizing the concepts of directional bias and
persistence as introduced by Weavers et al.^[228]33. Directional bias,
denoted by
[MATH: α :MATH]
, represents the angle between the instantaneous velocity of CTLs and
the direction from CTLs to the nearest point on the cancer-cell
cluster, reflecting the ability of cells to move in the direction of an
attractant source. Additionally, directional persistence, denoted by
[MATH: β :MATH]
, measures the angle between CTL velocities in two consecutive steps,
indicating the ability of cells to maintain their movement direction.
The bias angle distributions range from 0° to 180° whereas the
persistence angle distributions range from −180° to 180°. The cosine
function was applied to the angles to normalize the values between −1
and 1, where a higher value indicates a greater bias or persistence.
For each CTL spot, the distance to the nearest pixel of the cancer
cluster boundary was calculated. The trajectory segments of CTLs
located outside the cancer clusters were exclusively selected for
analyzing directional bias and persistence. These segments were then
grouped based on their distances from the nearest cancer-cluster
boundary, with each group defined by 5 μm intervals. Subsequently, the
averaged bias or persistence was computed for each distance-interval.
To compare the behavioral dynamics of CTLs, the average bias and
persistence distributions for CTLs cocultured with OVA cancer cells and
those with WT cancer cells were calculated over three specific
time-periods: 0–6 h, 6–12 h, and 12–18 h, respectively (Fig. [229]2F,
G). Additionally, directional bias and persistence analyses were
conducted per 3-h intervals, as shown in the supplementary information
(Supplementary Fig. [230]5B, C).
Mean square displacement of CTLs
To investigate the motility rules and compared the rules in different
time intervals, the mean square displacement was calculated to
characterize the diffusive behavior for CTLs. The calculation was
performed using the following equation^[231]59:
[MATH: MSDn△t=∑1N
−nl
i,i+n2N−n :MATH]
Here,
[MATH:
li,i+<
/mo>n :MATH]
represents the Euclidean distance of a CTL from frame
[MATH: i :MATH]
to frame
[MATH: i+n :MATH]
; N represents overall time periods for calculation; n represents the
time lag.
The MSD can also be described using a power-law form^[232]60:
[MATH: MSD(n△t
)=4D
△tα :MATH]
The trajectory segments of CTLs before reaching the cancer cell
clusters were extracted. Then, the logarithm of the MSD as a function
of the time was fitted, and the linear portions of the curve were
identified (Supplementary Fig. [233]5A). To ensure credibility and
accuracy of the calculations, only the first 10 frames of continuous
segments in the tracks were selected for the fitting process.
Quantification of collagen fibers
To investigate potential variations in collagen fiber distribution and
arrangement between WT and OVA cancer cell environments, collagen fiber
orientation were analyzed around cancer clusters in two co-culture
systems at 0 h and 6 h, respectively. For the analysis of fiber
orientation, CT-fire software^[234]61 was utilized to extract and
calculate the directional alignment angles between collagen fibers and
cancer boundaries (Supplementary Fig. [235]6A–C).
Correlation test among directional bias, directional persistence and step
size
To investigate whether there exists correlation among directional bias,
directional persistence and step size, Spearman correlation tests were
performed by per 5 μm within 0–6 h, 6–12 h, and 12–18 h, respectively
(Supplementary Fig. [236]7).
Agent-based models
To further investigate the impact of observed differences in T-cell
motility patterns on their searching efficiency and the accumulation on
cancer-cell clusters, we constructed 2D agent-based models that
incorporated four key features of T cells: dwell time, directional
bias, directional persistence and step size. In addition, computational
simulations offer the flexibility to manipulate the initial
distributions of CTLs and cancer clusters, which is challenging to
control under experimental conditions. To address this, we developed
two models: a single-cancer-cluster model (Supplementary Fig. [237]8)
and a multi-cancer-cluster model (Fig. [238]3), allowing us to explore
the potential effects of the initial distribution of CTLs and cancer
clusters.
The multiple-cancer-cluster model is based on the following
assumptions:
1. The model is implemented on a 2D square region with dimensions of
590 μm × 590 μm, matching the experimental scale. The model
includes CTLs and cancer-cell clusters. To simplify the modeling
process, the models disregard the mobility of cancer cells, as well
as the proliferation or death of both CTLs and cancer cells, and
the weak correlation among step size, directional bias and
persistence.
2. To simplify the model, the initial distribution of cancer-cell
clusters is identical to Series 1 in the experiments. Initially,
200 CTLs are randomly dispersed outside cancer-cell clusters within
the defined region. The simulation utilized parameters derived from
Series 1 in the WT and OVA experimental systems, including step
size, directional bias, directional persistence as a function of
distances towards cancer-cell clusters, and dwell time
distributions within 6–12 h, as significant differences were
observed in the directional bias, persistence, and dwell time
distribution of CTLs during this period (Fig. [239]2D–G). The step
size, directional bias and persistence distributions are updated
with a resolution of 5 μm distance from the nearest cancer-cell
boundary.
3. When T cells are outside of cancer-cell clusters, they undergo a
distance-dependent persistent or bias walk.
If we specifically focus on the influence of the differences in the
distribution of directional persistence, we can determine the
position of T cells at any given time by following the equations
below.
[MATH:
Txt
+1=Tx0+∑i=0
t+1dicos∑0iβi
:MATH]
[MATH:
Tyt
+1=Ty0+∑i=0
t+1d<
/mi>isin∑0iβi
:MATH]
[MATH:
Txt<
/mi> :MATH]
and
[MATH:
Tyt<
/mi> :MATH]
represents the coordinates of a CTL at time t.
[MATH: β0
:MATH]
represents the angle between the initial motion vector and the
positive x-axis, which is randomly set.
[MATH: dt
:MATH]
represents the step size from t to t + 1, sampled from the
step-size distributions. Here, we developed an algorithm to realize
the sampling process from a known distribution.
[MATH: βt
:MATH]
represents the turning angles from t to t + 1, sampled from
persistence distributions based on the distance towards cancer-cell
clusters.
Similarly, if we focus solely on the role of directional bias, we
can determine the position of CTLs using the following equations.
The motion of CTLs depends on searching for the nearest
cancer-cluster point at any given time, with the coordinates of
denoted as
[MATH: Cxt,C<
/mi>yt
:MATH]
.
[MATH: xt,yt=⟨C<
mrow>xt−Txt,Cyt−Tyt⟩<
/mo> :MATH]
[MATH:
αx<
mi>t=cos
−1xt∣
∣⟨xt,yt<
/msub>⟩∣∣ :MATH]
[MATH:
Txt
+1=Txt+dicos(αt<
/mrow>−α<
mi>xt
) :MATH]
[MATH:
Tyt
+1=Tyt+disin(αt<
/mrow>−α<
mi>xt
) :MATH]
[MATH:
αxt :MATH]
represents the angle between the motion vector and the positive
x-axis at time t.
[MATH: αt
:MATH]
represents bias angle from t to t + 1 sampled from bias
distributions with respect to the distance towards cancer-cell
clusters.
4. Once T cells reach cancer-cell clusters, they can remain for a
duration time t sampled from the dwell-time distributions. Each CTL
underwent 800 computational simulation steps.
5. Then, the effects of the four patterns within two systems were
explored individually by altering the directional bias distribution
or directional persistence or step-size distribution, while
maintaining consistency in the other two patterns (Table [240]1).
In terms of searching efficiency, we computed the mean arrival time
for all CTLs to initially reach cancer clusters in each simulation,
which is depicted as a data point in the violin plot generated
using the Violinplot function
([241]https://github.com/bastibe/Violinplot-Matlab) illustrated in
Fig. [242]3D.
6. In addition, we examined the number of T cells accumulating on
cancer-cell clusters over time to explore whether the motion
outside of cancer-cell or adhesion could affect the accumulation of
CTLs on cancer cell clusters. Each simulation was repeated 100
times.
Table 1.
Computational modeling of the effects of the four motility patterns
Simulation Step size Persistence Bias Dwell time
1 WT - - -
OVA - - -
2 - WT - -
- OVA - -
3 - - WT -
- - OVA -
4 - - - WT
- - - OVA
[243]Open in a new tab
Bulk RNA-seq data processing
The raw sequencing reads were trimmed to remove sequences with low
quality and short read lengths using Fastp (version 0.20.0).
Subsequently, HISAT2 was employed to map the high-quality reads onto
the mouse reference genome mm10. The read count matrix corresponding to
the genes and samples was generated using featureCounts. Next, the
count matrix was subjected to the R package DESeq2 to identify
differentially expressed genes with a log2(fold change) cutoff of 1 and
a false discovery rate (FDR) of 0.05 (Fig. [244]4A, E). Additionally,
we transformed the count matrix into a TPM matrix to compare the
special genes expression levels among different samples (Supplementary
Fig. [245]10A) and R package ggplot2 and pheatmap were used to
visualize (Fig. [246]4B, F). GO term enrichment analysis (MF, molecular
function, CC, cellular component, BP, biological process) and KEGG
pathway analysis were performed using the R package
clusterProfiler^[247]62 based on the up-regulated and down-regulated
genes, with significance determined by an FDR threshold of less than
0.05 (Fig. [248]4C, D, G, H). The batch effect between samples from the
Round 1 and Round 2 count matrices was corrected using R package sva.
Principal component analysis (PCA) was performed using the prcomp
function and visualize the first two principal components with the
ggplot2 package (Supplementary Fig. [249]9A). Gene set variation
analysis (GSVA) was conducted for specific gene set using the R package
GSVA (Supplementary Fig. [250]9C). Gene set enrichment analysis (GSEA)
was conducted for GO, KEGG and HALLMARK gene sets using the R package
GSEA (Supplementary Fig. [251]9H).
ScRNA-seq data processing
Raw gene expression matrices were generated for each sample using the
Cell Ranger (v6.1.2) pipeline coupled with the mouse reference genome
version mm10-2020-A. The output filtered gene expression matrices were
analyzed using R software (version 4.1.2) with the Seurat package
(version 4.1.1)^[252]63. A custom R script was used to combine the
expression data and metadata from all libraries corresponding to a
single batch. The expression data matrix was loaded into a Seurat
object along with the library metadata for downstream processing. The
percentage of mitochondrial transcripts for each cell (percent.mt) was
calculated and added as metadata to the Seurat object. Cells were
further filtered before dimensionality reduction (nFeature_RNA > 200 &
percent.mt<10%). Expression values were then scaled to 10,000
transcripts per cell and log-transformed (NormalizeData function,
normalization.method = “LogNormalize”, scale.factor = 10000). We
calculated variable features using the FindVariableFeatures function
with selection.method = “vst” and nfeatures = 2000. Effects of
variables (percent.mt) was estimated and regressed out using the
ScaleData function (vars.to.regress = “percent.mt”), and the scaled and
centered residuals were used for dimensionality reduction and
clustering.
Dimensionality reduction
To reduce the dimensionality of the dataset, the RunPCA function was
conducted with default parameters on linear-transformation scaled data
generated by the ScaleData function. Next, the ElbowPlot and DimHeatmap
functions were used to identify proper dimensions of the dataset (first
11 PCs).
Cell clustering and cluster identification
We first performed the FindNeighbors function before clustering, which
takes as input the previously defined dimensionality of the dataset
(first 11 PCs). Clustering was performed at varying resolution values
using the FindClusters function, and we chose a value of 0.7 for the
resolution parameter for the initial stage of clustering. After
non-linear dimensional reduction and projection of all cells into
two-dimensional space by UMAP, clusters were assigned preliminary
identities based on the expression of combinations of known marker
genes for major cell classes and types.
Sub-clustering and secondary clustering
We filtered out low quality clusters and mixed cell clusters (C1, C6,
C9, C16, C17). Sub-Clustering was performed on the two major cell types
(CD8 T Cells and Melanoma) data subset to resolve additional cell types
and subtypes. Secondary Clustering was performed using a standard
workflow similar to before but with different dimensionality (first 12
PCs) and a value of 0.5 for the resolution parameter. We filtered out
three clusters (C9, C10, C11), renamed all clusters from C1 to C9 (as
shown in Fig. [253]5A) and then performed other downstream analyses.
Differential expression genes (DEGs) identification and functional enrichment
Differential gene expression testing was performed using the
FindAllMarkers function in Seurat with the parameter only.pos = TRUE.
DEGs were filtered using a minimum log2(fold change) of 0.5, a minimum
expression percentage value of 0.3 and a maximal p value of 0.05. Dot
plots showing the relative average expression of top 10 DEGs in
subclusters were visualized (Fig. [254]5B). All DEGs’ relative average
expression in subclusters was visualized using the DoHeatmap function
in Seurat (Supplementary Fig. [255]12A). GO (BP) enrichment analysis
for the functions of the DEGs was conducted using the
clusterProfiler^[256]62 (version 4.2.2) R package (Supplementary
Fig. [257]12B–D).
DEGs analysis between two groups
A Wilcoxon rank-sum test for differential gene expression comparing
CD8T-c2 (C9) versus CD8T-cx (merged clusters of C5 to C8, Fig. [258]5C)
in OVActl was performed using the FindAllMarkers function
(min.cells.group = 0, min.pct = 0, logfc.threshold = -Inf) in Seurat
(Fig. [259]5D). Differential gene expression comparing
Melanoma_c3_Cxcl9 versus Melanoma_c4_Rhox5 (cell type of
secondary-clustered dataset) was also performed with the same function
and parameters (Fig. [260]5E). Percentage difference (Δ Percent of
Cells) and log2(fold change) are shown in a volcano plot. Genes
highlighted in red or blue have adjusted p < 0.05.
Cell-cell communication analysis (CellChat)
We used CellChat^[261]37 (v1.4.0) to infer cell-cell communication
based on prior ligand-receptor interaction databases. We separately
loaded the normalized counts data of WT group (WTcancer, WTctl) and OVA
group (OVAcancer, OVActl) extracted from Seurat object into CellChat
and followed the workflow recommended in CellChat to infer and
visualize cell-cell communication network. The database of known
receptor-ligand pairs was applied to assess cell-cell communication in
the selected target clusters (CD8T-c2, CD8T-cx, Melanoma_c3_Cxcl9,
Melanoma_c4_Rhox5) of OVA group (Fig. [262]5F and Supplementary
Fig. [263]12F). Interactions were trimmed based on significant sites
with p < 0.01.
Statistics and reproducibility
Statistical analysis of quantitative motility patterns and
computational modeling were conducted using MATLAB. Results are
presented as mean ± s.e.m., with individual data points shown, unless
otherwise stated. The statistical test methods and their corresponding
p-values are listed in Supplementary Data [264]8. For all relevant
figures, p-value annotations are as follows: ns: p ≥ 0.05; *:
0.01 < p < 0.05; **: 0.001 < p ≤ 0.01; ***: p ≤ 0.001, unless otherwise
indicated. To ensure reproducibility, at least five parallel
experiments were conducted, and statistical analyses were performed on
the resulting data (Fig. [265]2, Supplementary Figs. [266]3–[267]5).
Computational simulations were repeated 100 times. For clarity in
visualization, data points are shown with a sub-sampling interval of 5
(Fig. [268]3D–H, Supplementary Fig. [269]8D–G). To enhance
visualization clarity, the curves were smoothed using a smoothing
parameter of 0.05 (Figs. [270]2C, [271]3D–H, Supplementary
Figs. [272]3B, C, [273]8D–G). Violin plots display the distribution
quartiles, with the white dot indicating the median and individual
points indicating raw data (Fig. [274]3D, Supplementary Fig. [275]8C).
For bulk RNA-seq, samples were derived from multiple biological
replicates across at least two independent experiments. For single-cell
RNA-seq, pooled cells from ~40 wells per condition were used from a
single experiment. All group sizes (n) refer to biological replicates
and are indicated in the figure legends.
Supplementary information
[276]Supplementary Information^ (20.3MB, pdf)
[277]Supplementary Data 1^ (37.5KB, xlsx)
[278]Supplementary Data 2^ (2.5MB, xlsx)
[279]Supplementary Data 3^ (462.1KB, xlsx)
[280]Supplementary Data 4^ (5.2MB, xlsx)
[281]Supplementary Data 5^ (2.4MB, xlsx)
[282]Supplementary Data 6^ (13.1KB, xlsx)
[283]Supplementary Data 7^ (61KB, xlsx)
[284]Supplementary Data 8^ (19.1KB, xlsx)
[285]Supplementary Movie 1^ (8.5MB, avi)
[286]Supplementary Movie 2^ (8.5MB, avi)
[287]Supplementary Movie 3^ (7.5MB, mp4)
[288]Supplementary Movie 4^ (7.6MB, mp4)
[289]42003_2025_8568_MOESM14_ESM.pdf^ (118.8KB, pdf)
Description of Additional Supplementary Files
[290]nr-reporting-summary^ (2.6MB, pdf)
[291]Transparent Peer Review file^ (5.9MB, pdf)
Acknowledgements