Abstract
Meningioma is the most prevalent primary brain tumor with extensive
intra-tumoral heterogeneity and high recurrence rates, particularly in
high-grade meningiomas. Despite advancements in understanding the
molecular underpinnings of meningiomas, the longitudinal evolutionary
trajectory and cellular diversity of recurrent tumors remain elusive.
In this study, we perform single-nuclei sequencing of matched primary
and recurrent meningiomas to explore the dynamic transcriptional
heterogeneity and evolutionary trajectory of meningioma tumor cells, as
well as their molecular interactions with tumor-associated immune cells
that shape the complex milieu of the meningioma microenvironment. Our
findings reveal that both primary and recurrent meningiomas constitute
diverse cellular compositions and hierarchies, where recurrent tumor
cells are characterized by enrichments of cell cycle activities and
proliferative kinetics. Integrative RNA velocity and latent time
uncover divergent transcriptional trajectories in recurrent tumors,
demonstrating multidirectional transitions with the dominance of
COL6A3, which confers higher risk vulnerability and treatment
resistance. Tumor microenvironment analysis further reveals enrichments
of COL6A3-mediated interactions between immunosuppressive macrophages
and tumor cells in recurrent meningiomas. Collectively, these results
provide profound insights into the complex evolutionary process of
meningiomas and suggest potential therapeutic strategies for the
treatment of recurrent tumors.
Subject terms: Cancer genomics, CNS cancer, Cancer microenvironment
__________________________________________________________________
Tumour evolution and heterogeneity of recurrent meningioma tumours
remains to be explored. Here, the authors perform single nuclei
sequencing of matched primary and recurrent meningiomas and reveal
diverse cellular compositions and hierarchies.
Introduction
Meningioma is the most common primary central nervous system (CNS)
tumor, with a prevalence of 9.51 diagnoses per 100,000
individuals^[50]1. While the majority of meningiomas are classified as
grade I benign tumors, approximately 20–30% are diagnosed as either
atypical (grade II) or anaplastic (grade III) tumors, which exhibit
more aggressive and malignant phenotypes with high recurrence
rates^[51]2. This has led to a growing interest in investigating the
unique molecular properties of meningiomas to uncover the underlying
mechanisms of tumorigenesis and to identify potential therapeutic
targets that could lead to treatment innovation.
Previous studies employing integrative genomic and transcriptomic
approaches have provided profound insights into the molecular
underpinnings of meningioma biology^[52]3–[53]14. These efforts have
resulted in the development of a new molecular classification system
that categorizes meningiomas into distinct subtypes, each characterized
by unique genomic alterations. These subgroups comprise immunogenic
(MG1), NF2 wild-type (MG2), hypermetabolic (MG3), and proliferative
(MG4), each defined by the activation of distinct biological pathways,
such as immune response, angiogenesis, lipid metabolism, and cell cycle
regulation.
In parallel, additional studies have demonstrated that meningiomas
exhibit a high degree of plasticity and undergo subtype transitions
between multiple cellular states^[54]15,[55]16. However, most existing
analyses have predominantly focused on primary meningiomas, leaving the
complex cellular hierarchy and dynamics of recurrent meningiomas
relatively unexplored and therapeutically unresolved due to the limited
understanding of their evolutionary trajectory. The molecular and
transcriptional heterogeneity within the malignant cell populations
provides meningiomas with diverse functional states that enable them to
evade treatment intervention, leading to tumor recurrence and
progression.
To address these unresolved aspects, we perform single-nuclei RNA
sequencing (snRNA-seq) of paired primary and recurrent tumors from
meningioma patients. We construct the longitudinal evolutionary
trajectory of meningiomas to examine the tumor cellular states over
time, identifying risk factors and potential therapeutic targets
driving the malignant state of recurrent meningiomas. In addition, by
analyzing the cellular composition of TME, including macrophages and
lymphocytes, we provide comprehensive insights into the dynamic
cellular interactions that reprogram the transcriptional states of
malignant tumor cells.
Results
Longitudinal profiling of meningioma at single-cell resolution
To explore the longitudinal evolutionary trajectory of meningioma, we
collected 14 matched fresh-frozen primary and recurrence samples from
seven meningioma patients (Fig. [56]1a). Twelve matched pairs were
obtained from primary and first recurrent tumors, while one pair
consisted of the first and second recurrences. Three patients were
diagnosed with grade I meningothelial or fibrous meningiomas, and the
remaining four were diagnosed with either grade II atypical or grade
III anaplastic meningiomas. All patients underwent craniotomy for tumor
removal, with an average age at diagnosis of 50 years. All 14 tumor
specimens were profiled using a droplet-based 10x snRNA-seq platform,
resulting in a total of 100,175 cells. Of these, 74,979 cells passed
the stringent quality control, with 36,601 genes detected in total
(Supplementary Table S[57]1, [58]2). After batch correction, cells were
assigned into 11 distinct clusters, including tumor cells (37,460
cells), lymphocytes (10,779 cells), myeloid cells (24,964 cells),
endothelial cells (714 cells), and additional minor populations based
on high variability in gene expression profiles (Fig. [59]1b–d and
Supplementary Table S[60]3, [61]4). Tumor cells were further
distinguished and verified based on copy-number variation (CNV)
analysis (Fig. [62]1e and Supplementary Fig. [63]S1). Tumor cells
exhibited a high degree of diversity at both intra- and inter-tumoral
levels, underscoring the heterogeneous nature of meningioma tumor
cells^[64]17,[65]18. Non-tumor cells were predominantly intermixed
across patients without notable batch effects. Tumor cells comprised
the majority of the cell populations, accounting for 53.3% (19,112 out
of 35,881) and 46.9% (18,348 out of 39,098) in primary and recurrent
tumors, respectively (Fig. [66]1f). The next largest populations were
macrophages (21,473 cells), followed by lymphocytes (10,479 cells). The
distribution of cell types varied among patients, with dendritic cells
predominantly found in primary tumors, whereas macrophages and
lymphocytes were more prevalent in recurrent tumors.
Fig. 1. A single-cell atlas of longitudinal meningioma.
[67]Fig. 1
[68]Open in a new tab
a Workflow for sample collection, sequencing and analysis of 14
meningioma. Created in BioRender. Lab, S. (2025)
[69]https://BioRender.com/pzifbd8b tSNE visualization of all 74,979
cells from 14 meningioma samples with cells colored according to the
corresponding cell types. c Heatmap for expression of marker genes in
all cell types of meningioma (left) and tSNE visualization of each
marker genes (right). d tSNE visualization from 14 meningioma with
cells colored according to the group, grade, patient, and sample. e CNV
plot that classifies tumor cells using immune cells as a reference
through inferCNV (MEN01-R sample). f Cell type proportion of meningioma
samples on a per-sample basis. Source data are provided as a Source
Data file.
Longitudinal evolutionary trajectory reveals variability in RNA velocity and
increased proliferative states in recurrent meningiomas
To identify biological features constituting the unique transcriptional
states of recurrent meningiomas, we performed genome-wide differential
analysis between primary and recurrent tumor cells (Fig. [70]2a). As a
result, we identified several transcriptomes actively enriched in
recurrent meningiomas, including DNA repair encoding molecules (POLQ,
BRIP1), cell cycle regulators (FOXM1) and extracellular matrix-related
genes (COL6A3) (Fig. [71]2b). Pathway enrichment analysis consistently
indicated that recurrent tumors were characterized by the activation of
genes functionally involved in cell cycle activity, proliferation
kinetics, DNA repair mechanisms, and metastasis (Fig. [72]2c). In
contrast, primary tumors exhibited increased expression of APOE, SOD3,
and HSPA6, which were involved in hypoxia, metabolism activity, and TNF
signaling. CNV analysis highlighted consistent loss in ARID1A and NF2
genes in both primary and recurrent meningiomas, while amplification of
angiogenesis-related genes such as ROBO1 and ROBO2 was more prevalent
in recurrent tumor cells (Fig. [73]2d, Supplementary Fig. [74]S2A and
Supplementary Table S[75]5). RNA velocity analysis using
velocyto^[76]19 revealed distinct gene expression dynamics between
primary and recurrent meningiomas. Recurrent tumors showed significant
variability in RNA velocity and latent time, characterized by
multidirectional transitions and increased expression of B2M in the
early stage, followed by a replacement with COL6A3 in the later stage
(Fig. [77]2e). Conversely, primary meningiomas demonstrated relatively
stable velocity with a unidirectional transition from increased
expression of CCND2 to LRP1B.
Fig. 2. Transcriptional heterogeneity and trajectory of primary and recurrent
meningioma tumor cells.
[78]Fig. 2
[79]Open in a new tab
a tSNE visualization of all 37,460 tumor cells from 14 meningioma
samples colored according to the primary or recurrence status. b
Differential gene expression analysis between primary and recurrent
tumor cells (19,112 primary tumor cells and 18,348 recurrent tumor
cells from 14 meningioma samples). p-values are from the two-sided Wald
test. Source data are provided as a Source Data file. c Pathway
enrichment analysis based on differentially expressed genes in primary
and recurrent tumor cells (19,112 primary tumor cells and 18,348
recurrent tumor cells from 14 meningioma samples). p-values are from
two-sided Fisher’s exact test. d Copy number alteration analysis of
primary and recurrent tumor cells from the inferCNV results (19,112
primary tumor cells and 18,348 recurrent tumor cells from 14 meningioma
samples). e RNA velocity, latent time and RNA phase portraits
(unspliced versus spliced) colored by cluster in primary tumor cells
(n = 19,112) (top) and recurrent tumor cells (n = 18,348) (bottom). f
Cell type proportion of tumor cells based on molecular subtypes between
primary (n = 19,112) and recurrent (n = 18,348) tumor cells. Detailed
cell counts are available in the Source Data file. g Cell type
proportion of cell cycle states based on molecular subtypes between
primary (n = 19,112) and recurrent (n = 18,348) tumor cells. Source
data are provided as a Source Data file. h Pathway enrichment score
across four molecular subtypes in tumor cells (19,112 primary tumor
cells and 18,348 recurrent tumor cells from 14 meningioma samples). The
colors represent the statistical significance (p-value), and the size
of the node represents the number of genes in the corresponding
pathways. i, j Two-dimensional representation of cellular subgroup of
tumor cells (n = 37,460). Yellow dots represent primary tumor cells
(n = 19,112), and green dots represent recurrent tumor cells
(n = 18,348). Each quadrant corresponds to one of the molecular
subtypes. The position of tumor cells (dots) reflects their relative
meta-modules scores, which were derived from ssGSEA analysis of the
four molecular subtypes and their colors reflect the latent time. Pie
chart show proportion of latent time group (early, intermediate, and
late) in each subgroup. Source data are provided as a Source Data file.
Previous studies have defined four clinically applicable molecular
classifications of meningiomas, each associated with distinct clinical
outcomes and unique genomic aberrations. MG1 tumors demonstrated the
most favorable prognosis with activation of immunogenic
characteristics. MG2 tumors were marked by the absence of NF2 genomic
alterations, while MG3 and MG4 subtypes were defined by increased
hypermetabolic and proliferative features, respectively^[80]10. We
adopted these subtype transcriptional signatures to annotate all tumor
cells in primary and recurrent meningiomas (Supplementary Table
S[81]6). Among the four identified molecular subtypes, MG2 tumor cells
were characterized as NF2 wild-type in both primary and recurrent
tumors, as determined by inferCNV analysis (primary CNV score: 0.15,
recurrent CNV score: 0.05). In contrast, the remaining subtypes
collectively exhibited NF2 loss, reinforcing that the
transcriptome-based subtype classification aligns with CNV alterations
(Supplementary Fig. [82]S2B). Notably, primary tumors demonstrated a
predominance of the MG1 immunogenic subtype, whereas recurrent tumors
exhibited a higher proportion of MG2 (NF2 wild-type) tumor cells
(Fig. [83]2f). In addition, cell cycle analysis revealed a greater
prevalence of proliferating cells in the S and G2M phases among
recurrent tumor cells compared to primary tumors across all four
subtypes (chi-square test P < 0.001, Fig. [84]2g). To investigate
subtype transitions upon recurrence, we assigned the most predominant
subtype within each tumor as its representative classification.
Interestingly, the patient experiencing a second recurrence exhibited a
pronounced shift from MG1 and MG2 to the MG3 hypermetabolism subtype,
which was significantly associated with poor clinical outcomes
(Supplementary Fig. [85]S3B). Functional analysis of each molecular
subtype demonstrated similarities in pathway activities between primary
and recurrent tumors. MG1 subtypes showed elevated levels of various
immunogenic signaling pathways, while MG2 tumors exhibited increased
angiogenesis-related and hypoxia pathways (Fig. [86]2h). The MG3
subtype was characterized by the upregulation of metabolic pathways
such as adipogenesis and fatty acid metabolism. In contrast, MG4 cells
showed enrichment in cell cycle and FOXM1-related pathways, which have
been previously recognized as major drivers in high-grade
meningiomas^[87]20,[88]21. Notably, differences in pathway activities
between primary and recurrent meningiomas were observed for the MG3 and
MG4 subtypes. MG3 recurrent tumors also showed activation of the cell
cycle and mismatch repair activities, while MG4 primary tumors
exhibited higher expression of the TCA cycle and adipogenesis. When we
compared the subtype composition based on WHO grades in recurrent tumor
cells, we observed a higher proportion of immunogenic tumor cells in
low-grade tumors, while high-grade meningiomas were marked by
enrichments of hypermetabolic and proliferative tumor cells
(Supplementary Fig. [89]S4A, B). Differentially expressed gene analysis
highlighted activation of FOXM1 in high-grade meningiomas, aligning
with previous findings (Supplementary Fig. [90]S4C). Pathway analysis
consistently demonstrated enrichments of cell-cycle and EMT-related
pathways in high-grade tumors. In contrast, low-grade meningiomas
showed increased immune-related signals such as inflammatory response
and IL6 pathways (Supplementary Fig. [91]S4D). Cell cycle analysis
further revealed a higher proportion of tumor cells in S and G2M phases
in high-grade meningiomas compared to low-grade tumors (Supplementary
Fig. [92]S3E). We also conducted a similar analysis on primary tumor
cells, observing largely consistent pathway enrichment trends between
primary and recurrent tumors (Supplementary Fig. [93]S4F–H). However,
we observed few key differences, particularly in pathways related to
EMT and DNA repair, which were significantly enriched in the recurrent
high-grade tumors but not in the primary high-grade meningiomas.
Integration of RNA velocity and latent time with subtype
transcriptional signatures revealed that early-state tumor cells were
predominantly observed in the immunogenic subtype (MG1) in both primary
and recurrent meningiomas, while late-state cells were more prevalent
in hypermetabolic (MG3) and proliferative (MG4) subtypes for primary
and recurrent tumors, respectively (Fig. [94]2i, j and Supplementary
Fig. S[95]3A). Collectively, these findings revealed distinct cellular
trajectories between primary and recurrent meningiomas and highlighted
differences in molecular subtype composition that define the unique
cellular profiles over the course of the disease.
Identification of COL6A3 as a major driver of risk vulnerability and
potential therapeutic target in recurrent meningioma
Previous studies have developed risk-associated transcriptional
signatures that could aid in predicting clinical outcomes in meningioma
patients using a large multi-institutional cohort^[96]22. These
polygenic biomarkers successfully discriminated between patients’
outcomes and demonstrated enhanced performance in predicting treatment
response to radiotherapy. Leveraging both “high-risk” and “low-risk”
gene signatures from prior studies, we assessed the risk states of
primary and recurrent meningioma tumor cells. Notably, the “high-risk”
signature was significantly upregulated in recurrent tumors, whereas
the “low-risk” signature was more prevalent in primary tumors
(Fig. [97]3a). Furthermore, the risk transcriptional signatures
significantly correlated with the latent time model, where the
“high-risk” score showed a strong correlation and the “low-risk”
inversely correlated (Supplementary Fig. [98]S3D). This association was
particularly pronounced in primary tumor cells, where late-latent tumor
cells exhibited elevated “high-risk” activity, while early-latent tumor
cells were predominantly associated with “low-risk” signature
(Fig. [99]3b). To identify molecular determinants driving increased
clinical risk in recurrent meningiomas, we developed an elastic-net
regression model, incorporating transcriptome data with latent time and
risk signature activity. As a result, we identified a set of candidate
targets, including COL6A3, NPW, and NRGN that were significantly
associated with higher risk indications (Fig. [100]3c). Among these,
COL6A3 (collagen type VI alpha 3 chain), previously recognized as an
essential component of extracellular matrix (ECM)^[101]23–[102]26,
emerged as one of the most robust hits. Prior studies have shown that
COL6A3 promoted increased resistance to platinum-based treatments, such
as cisplatin and oxaliplatin, in ovarian and breast
tumors^[103]27–[104]29. Moreover, ColVI chain-encoding molecules were
highly expressed in meningiomas^[105]30. Comparison of COL6A3 mRNA
expression levels based on risk signature activities revealed that
tumor cells with high “high-risk” scores demonstrated increased COL6A3
expression, whereas the opposite trend was observed for “low-risk”
signature scores (Fig. [106]3d). Furthermore, COL6A3 exhibited a
significant correlation with latent time, particularly in recurrent
meningioma tumor cells (Supplementary Fig. [107]S3E). When categorized
into four distinct molecular subtypes of meningiomas (MG1, MG2, MG3,
and MG4), we discovered that tumor cells with high COL6A3 expression
levels were primarily observed in recurrent tumor cells, specifically
in MG4 (proliferative) subgroup, followed by MG3 (hypermetabolism) and
MG2 (NF2 wild-type) (Fig. [108]3e).
Fig. 3. Risk signature analysis of primary and recurrent tumor cells
identified COL6A3 as a driver of risk vulnerability.
[109]Fig. 3
[110]Open in a new tab
a Risk score (“high-risk” and “low-risk”) comparison between primary
and recurrent meningioma tumor cells (n = 37,460). Box plots show the
median, interquartile range, and 1.5 × IQR whiskers; outliers are shown
as individual points. Statistical significance was assessed using the
two-sided Wilcoxon rank-sum test. Source data are provided as a Source
Data file. b Comparison of higher risk score (left) and lower risk
score (right) in primary tumor cells categorized by early and late
latent time groups (n = 5,733). Box plots show the median,
interquartile range, and 1.5 × IQR whiskers; outliers are shown as
individual points. Statistical significance was assessed using the
two-sided Wilcoxon rank-sum test. Source data are provided as a Source
Data file. c Top genes associated with latent time identified through
the elastic-net regression model. The x-axis indicates the correlation
between gene expression and latent time across individual cells, and
the y-axis shows the frequency of selection across 100 bootstrap
iterations of the elastic-net model. p-values were calculated by
Spearman correlation tests. Source data are provided as a Source Data
file. d Expression level of COL6A3 gene according to risk score groups.
Each group of tumor cells were selected based on “high-risk” (top) or
“low-risk” (bottom) signatures and termed “top 25% tumors” and “bottom
25% tumors” based on their risk signature activities. e Log[2]
normalized gene expression level of COL6A3 between primary and
recurrent tumor cells based on molecular subtype (two-sided Wilcoxon
rank-sum test). Asterisks indicate statistical significance: p < 0.0001
(****). f Violin plot showing the comparison of COL6A3 expression
levels between primary and recurrent samples in the validation cohort
(n = 112 samples). Box plots show the median, interquartile range, and
1.5× IQR whiskers; outliers are shown as individual points. Statistical
significance was assessed using the two-sided Wilcoxon rank-sum test. g
RFS (Relapse-Free Survival) according to COL6A3 expression level in the
validation cohort (n = 110). The “high” group represents patients with
COL6A3 expression that’s higher than the threshold determined by the
maxstat algorithm, and the “low” group depicts patients with COL6A3
expression lower than the threshold. Source data are provided as a
Source Data file. h Bar plot comparing cell proliferation between
siRNA-COL6A3 and control in IOMM-Lee and SF3061 meningioma cell lines.
Each bar represents the mean of n = 4 biological replicates, and the
error bars indicate the standard deviation (SD). Source data are
provided as a Source Data file. i Bar plot comparing the cell cycle
proportion between siRNA-COL6A3 and control in IOMM-Lee and SF3061 cell
lines. Each bar represents the mean of n = 6 biological replicates, and
the error bars indicate the standard deviation (SD). Source data are
provided as a Source Data file. j IHC staining images of COL6A3 in
matched primary (MEN08, MEN09) and recurrent (MEN08, MEN09)
meningiomas. Scale bar, 300 μm. k Violin plot showing the paired
comparison of COL6A3 protein expression levels between matched primary
and recurrent samples through IHC staining. Box plots show the median,
interquartile range, and 1.5 × IQR whiskers; outliers are shown as
individual points. Statistical significance was assessed using the
two-sided Wilcoxon rank-sum test. A total of 18 primary and 18
recurrent samples were analyzed, including 12 matched pairs. Source
data are provided as a Source Data file.
Moreover, we analyzed an independent external cohort consisting of 112
meningioma RNA-seq samples from both primary and recurrent
tumors^[111]14. This analysis revealed a significantly higher
expression of COL6A3 in recurrent samples compared to primary tumors
(Wilcoxon test, p = 0.0054) (Fig.[112]3f). To further assess whether
COL6A3 contributes to increased recurrence vulnerability in
meningiomas, we analyzed another external cohort of 110 meningioma
patients with both RNA sequencing data and clinical outcomes^[113]31.
Notably, patients with high COL6A3 mRNA expression exhibited
significantly worse relapse-free survival (RFS) probabilities
(Fig. [114]3g). Functional validation through in vitro assays further
supported these findings, as siRNA-mediated silencing of COL6A3
significantly reduced cellular proliferation in IOMM-Lee and SF3061
meningioma cell lines (Fig. [115]3h). In addition, COL6A3 knockdown led
to a marked reduction in cell cycle activity, highlighting its
essential role in tumor progression and proliferation (Fig. [116]3i).
To corroborate these results, we performed immunohistochemical (IHC)
staining of COL6A3 in 18 primary and 18 recurrent samples (including 12
matched patient pairs). Consistent with our prior findings, recurrent
tumors exhibited significantly higher expression levels of COL6A3
compared to primary tumors (Fig. [117]3j, k). Furthermore,
transcriptome analysis across three independent external cohorts
([118]GSE183653, [119]GSE197763, and [120]GSE136661) confirmed that
COL6A3 expression was significantly elevated in recurrent tumors
compared to primary, as well as in high-grade meningiomas
(Supplementary Fig. [121]S5A, B)^[122]14,[123]32. Collectively, these
findings underscore COL6A3’s pivotal role in tumor progression and
recurrence. Given its contribution to tumor cell proliferation and
recurrence, COL6A3 represents a promising therapeutic target for
high-risk recurrent meningiomas.
TME analysis reveals cell-cell interaction between tumor cells and C1Q +
macrophages that drive COL6A3 crosstalk in recurrent meningiomas
Next, we investigated whether cell-cell interactions between tumor
cells and immune cells in the tumor microenvironment (TME) influence
the cellular diversity between primary and recurrent tumors. Focusing
on myeloid cells, the primary source immune cell population in the
brain, we delineated the meningioma surrounding microenvironment.
Myeloid cells were classified into eight distinct clusters based on
unique gene expression profiles (Fig. [124]4a). Utilizing previously
established transcriptional signatures^[125]10,[126]15,[127]33–[128]35,
myeloid cell subpopulations were further identified as C1Q+ macrophages
(C1QA, C1QB, AIF1), M1 macrophages (STAT1, NFKB1, CXCL10), M2
macrophages (CD163, MRC1), proliferative macrophages (BRIP1, MKI67),
and other defined populations. (Fig. [129]4b and Supplementary
Table.S[130]7). Among them, we discovered a high prevalence of
inflammatory macrophages and monocytes in primary tumors, whereas
recurrent tumors exhibited higher infiltration of C1Q + , M1, M2, and
proliferative macrophages (P < 0.001; Fig. [131]4c). Notably, the
infiltration of C1Q + tumor-associated macrophages (TAMs) has been
associated with poor prognosis, promoting T cell exhaustion and immune
tolerance induction in healthy tissues^[132]36–[133]38. Pathway
analysis revealed distinct features between primary and recurrent
myeloid cells, where hemostasis regulation, ECM organization, and cell
cycle pathways were prominently activated in recurrent tumors
(Fig. [134]4d). Primary tumors, on the other hand, showed activation of
pathways related to cell migration and anti-tumor responses, such as
lymphocyte activation and tumor necrosis factors. To account for
potential confounding sources of variation between primary and
recurrent tumors, we adopted differential cell abundance analysis using
k-nearest neighbor (KNN) based graphs. Specifically, we utilized miloR
to identify transcriptionally similar cell neighborhoods within
macrophage and monocyte populations. These results further supported
our findings, showing enrichments of inflammatory macrophages and
monocytes in primary tumors, while C1Q +, M1, and M2 macrophages were
identified as the primary source of TAM infiltration in recurrent
tumors (Fig. [135]4e). These findings highlight that recurrent tumors
are characterized by a more pronounced pro-tumorigenic potential
through the infiltration of M2 and C1Q + macrophages.
Fig. 4. TME analysis reveals cell-cell interaction between tumor cells and
myeloid cells in longitudinal meningioma.
[136]Fig. 4
[137]Open in a new tab
a tSNE visualization of 21,800 cells from macrophages and monocytes
colored according to corresponding cell types. Detailed cell counts are
available in the Source Data file. b Heatmap for expression of
representative marker genes in myeloid cell types of meningioma
(n = 21,800 cells)(top) and tSNE visualization of each marker genes
(bottom). c Cell type proportion of myeloid cell (n = 21,800) by
per-sample basis. Detailed cell counts are available in the Source Data
file. d The pathway enrichment analysis between primary (n = 10,776)
and recurrent (n = 11,024) myeloid cells. p-values are from two-sided
Fisher’s exact test. e Beeswarm plot showing the distribution of log[2]
fold change in abundance between conditions in neighborhoods from
different cell type clusters (n = 21,341 cells). Each dot represents a
neighborhood, a group of cells clustered based on shared gene
expression profiles. f Comparison of integral signal with input and
output signals between myeloid (n = 21,800) and tumor cells
(n = 37,460). Statistical significance was assessed by permutation in
CellChat. g Circle plots showing the interactions and strengths among
different cell types in primary (left) (n = 10,776) and recurrent
(n = 11,024) myeloid cells (right). h IHC staining images of COL6A3,
C1QA, and CD163 in matched primary (MEN10) and recurrent (MEN10)
meningiomas. Representative images are shown; IHC performed on 36 tumor
samples (18 primary meningioma and 18 recurrent meningioma samples).
Scale bar, 300 μm. i Scatter plot showing correlation between COL6A3
and C1QA expression scores in primary (n = 18) and recurrent (n = 18)
meningiomas. p-values were calculated by Spearman correlation tests.
Source data are provided as a Source Data file.
We speculated that the observed diversity in tumor cell populations
between primary and recurrent tumors might result from cell-cell
interactions within the TME, which could reprogram the transcriptional
state of tumor cells. To explore these associations, we performed
ligand-receptor analysis using CellChat to quantitatively assess the
cellular crosstalk between tumor cells and myeloid cells. In primary
tumors, the BMP5-BMPR1a signaling axis was significantly enriched
between tumor cells and myeloid cells across various subtypes
(Fig. [138]4f, g). Given its prominent role in prostate cancer, where
BMPR1a promotes tumor growth while its loss impairs tumor progression,
primary meningiomas are expected to exhibit similar pro-tumorigenic
characteristics mediated by BMPR1a signaling. In addition, BMP
signaling inhibition has been shown to disrupt M2 macrophage
polarization, potentially supporting an anti-tumorigenic tumor
microenvironment^[139]39,[140]40. In recurrent tumors, significant
enrichment of COL6A3-CD44 and COL6A3-ITGA9/ITGB1 crosstalk between
tumor cells and myeloid cells were observed. Particularly, C1Q+
macrophages were enriched in the MG2, MG3, and MG4 subtypes,
corroborating our prior findings that COL6A3 activation drives
increased risk vulnerability in recurrent meningiomas. Interestingly,
immunogenic recurrent meningiomas did not exhibit any significant
crosstalk in the COL6A3 signaling axis (Fig. [141]4f, g). To further
investigate the role of COL6A3 with C1Q + macrophages in shaping the
recurrent meningioma microenvironment, we performed IHC staining on
primary and recurrent meningiomas to assess the spatial distribution of
COL6A3, C1QA (a marker of C1Q + macrophages), and CD163 (a marker
associated with immunosuppressive macrophages) (Fig.[142]4h, i and
Supplementary Fig. [143]S5C, D). IHC analyses revealed significantly
higher expression levels of COL6A3 and C1QA in recurrent tumors
compared to primary tumors at the protein level. This difference was
particularly pronounced in matched primary-recurrent pairs, further
supporting its role in tumor recurrence (Supplementary Fig. [144]S5D).
Moreover, at the individual tumor level, COL6A3 expression was
significantly correlated with C1QA, suggesting a potential interaction
between COL6A3-expressing tumor cells and C1Q+ macrophages in the
recurrent tumor microenvironment. Repolarization of M1 into M2
macrophages is mediated by multiple integrin signals, particularly
integrin β7, integrin α9, and integrin β1 in a TGFβI-dependent
manner^[145]41. We speculate that paracrine crosstalk between tumor
cells and C1Q + macrophages through the COL6A3 signaling axis promotes
conversion of tumor cells to a more malignant transcriptional state in
recurrent meningiomas. Taken together, these findings highlight the
critical role of COL6A3-mediated signaling in tumor progression and
immune modulation, underscoring its potential as a therapeutic target
for high-risk recurrent meningiomas.
Cellular architecture of lymphocytes in longitudinal meningioma
In addition to myeloid cells, we further assessed the cellular
architecture of lymphocytes. We identified 14 distinct lymphocyte
populations, including Natural Killer (NK) cells, CD4, and CD8 T cells,
which were further subclassified based on their functional states such
as CD4 regulatory T cells (CTLA4, FOXP3, ICOS), cytotoxic CD8 T cells
(CCL5, CCL4, GZMA), NKG2A + NK cells (KLRC1, GNLY,
KLRF1)^[146]42–[147]44 (Fig. [148]5a, b and Supplementary Table
S[149]8). Interestingly, recurrent meningiomas were characterized by
increased populations of CD4 naïve and NKG2A + CD8 T cells, while CD4
regulatory T cells and HSP + CD4 T cells were dominant in primary
tumors (Fig. [150]5c). Notably, infiltration of NKG2A + NK cells was
evident in high-grade recurrent meningiomas, while low-grade tumors
showed increased levels of GZMK + CD8 T cells in both primary and
recurrent meningiomas. Differential gene expression analysis further
revealed increased expressions of GZMB and S100A8 in recurrent tumors,
while heat shock protein-related genes, including HSPA6, HSPA1A, and
DUSP1, were more abundant in primary tumors (Fig. [151]5d).
Fig. 5. Distinctive characteristics of T cells and NK cells between primary
and recurrent meningiomas.
[152]Fig. 5
[153]Open in a new tab
a tSNE visualization of 10,491cells from T cells and NK cells colored
according to corresponding cell types. Detailed cell counts are
available in the Source Data file. b Dot plot for expression of marker
genes in T cell and NK cell types of meningioma (n = 10,491cells) (top)
and tSNE visualization of each marker gene (bottom). c Cell type
proportion of T cell (n = 9395) and NK cell (n = 1096) by patient and
grade. Detailed cell counts are available in the Source Data file. d
Differential gene expression comparing cells from NK/T cells between
primary and recurrent meningiomas (3434 primary cells and 7057
recurrent cells). P-values are from the two-sided Wald test. Source
data are provided as a Source Data file. e The pathway enrichment score
of NK/T cells from primary (n = 3434) and recurrent (n = 7057) tumors.
p-values are from two-sided Fisher’s exact test. f Violin plots of
ssGSEA-based pathway scores across NK/T cells (n = 10,491). Box plots
show the median, interquartile range, and 1.5 × IQR whiskers; outliers
are shown as individual points. Statistical significance was assessed
using the two-sided Wilcoxon rank-sum test. Source data are provided as
a Source Data file. g Comparison of Integral signal with input and
output signals between NK/T cells (n = 10,491) and recurrent MG1 tumor
cells (n = 3176) (left) and chord diagram showing the network of HLA
signaling pathways mediated only by HLA-E – CD94:NKG2A in different
cell types(right). Dot color reflects communication probabilities, and
dot size represents computed p-values. Statistical significance was
assessed by permutation test implemented in CellChat.
Gene Set Enrichment Analysis (GSEA) highlighted the activation of
metabolism-related pathways, including glucose deprivation in
lymphocytes from recurrent tumors, whereas primary tumors exhibited
enrichments of immune memory-related pathways (Fig. [154]5e). In
addition, lymphocytes from primary meningiomas exhibited adaptation to
a hypoxic environment, while recurrent tumors showed reduced lymphocyte
progenitors and diminished immune characteristics (Fig. [155]5f).
Cell-cell interaction between tumor cells and lymphocytes demonstrated
significant enrichments of crosstalk only within MG1 (immunogenic)
recurrent meningiomas (Fig. [156]5g). Notably, immunogenic tumor cells
interacted with various T cell populations through the SPP1-CD44 and
SPP1-ITGA4/ITGB1 signaling axis. SPP1 has been identified as a
prognostic marker in intrahepatic cholangiocarcinoma (ICC) and is
associated with poor survival^[157]45. We speculated that recurrent
tumor cells may inhibit T cell proliferation through SPP1-CD44
interactions, contributing to tumor progression. In addition,
activation of the HLA-E pathway was evident in the cell-cell
interaction between tumor cells with cytotoxic CD8 T cells, NKG2A + CD8
T cells, and NKG2A + NK cells. NKG2A has been regarded as an essential
immune checkpoint molecule, and prior studies have shown that tumor
cells elicit CD94 and NKG2A expressions on CD8 + T cells, leading to T
cell exhaustion^[158]46. In addition, inhibition of HLA-E:CD94-NKG2A
prevented pancreatic tumor metastasis through blood
circulation^[159]47. In recurrent meningiomas, the increased presence
of NKG2A + T cells suggests a unique tumor microenvironment that
potentially enhances tumor malignancy and risk vulnerability through
the HLA-E-NKG2A signaling axis. Collectively, these findings highlight
the diverse roles of lymphocytes in shaping the unique tumor
microenvironments of longitudinal meningiomas.
Discussion
Tumor recurrence has been a long-standing challenge in the treatment of
malignant brain tumors, including meningiomas^[160]14,[161]48. While
the majority of meningioma diagnoses are grade I tumors, over 20% of
the patients are diagnosed with high-grade meningiomas and exposed to
treatment resistance and increased vulnerability to tumor relapse.
Particularly, the longitudinal evolutionary trajectory of meningioma
remains poorly understood, especially at single-cell resolution.
In this study, we constructed a single-cell atlas of longitudinal
meningiomas to uncover the dynamic cellular composition and
transcriptional heterogeneity of meningioma tumor cells and
tumor-associated immune cells undergoing tumor evolution. We
investigated the distinct transcriptional trajectory between primary
and recurrent tumor cells, identified a driver of risk vulnerability
and potential therapeutic target, and explored cell-cell interactions
between tumor cells and tumor-associated immune cells within the tumor
microenvironment. Our findings revealed that recurrent meningiomas
exhibited activation of cell cycle and DNA repair pathways compared to
matched primary tumors. Moreover, RNA velocity and latent time models
indicated a predominance of hypermetabolic and proliferative tumor
cells during late phases in recurrent tumors, driven by increased
expression of B2M and COL6A3 during early and late phases,
respectively. By assessing the distinct longitudinal trajectory between
primary and recurrent tumors, we developed an elastic-net regression
model to identify potential therapeutic targets driving the risk
vulnerability of recurrence tumors. Notably, COL6A3 emerged as one of
the most robust hits, which we further validated using external RNA-Seq
cohorts of meningiomas patients and functional experiments. These
findings were consistent with previous reports in colorectal cancer,
pancreatic ductal adenocarcinoma, and glioblastoma, where COL6A3 has
been implicated in tumor progression, extracellular matrix remodeling,
and immune evasion, further supporting its pro-tumorigenic role across
multiple cancer types^[162]49–[163]51. In addition, by exploring
cellular architecture in the TME, we discovered dynamic cell-cell
interactions between tumor cells and immune cells, including
macrophages and T cells that shape the unique tumor surroundings of
primary and recurrent meningiomas. Specifically, recurrent meningiomas
were characterized by infiltration of immunosuppressive macrophages,
including C1Q+ and M2 macrophages. Ligand-receptor analysis revealed
that COL6A3-CD44 and COL6A3-ITGA9/ITGB1 crosstalk between tumor cells
and C1Q+ and M2 macrophages were highly enriched in recurrent tumors.
Functional validation through immunohistochemistry (IHC) confirmed a
significant upregulation of COL6A3 and C1QA in recurrent tumors, with a
strong correlation between COL6A3-expressing tumor cells and C1Q+
macrophages. These findings suggest that COL6A3 interacts with C1Q+
macrophages via ITGA9 + ITGB1 integrins, shaping the distinct immune
landscape of recurrent meningiomas. Integrin α9β1 has been known to
mediate interactions with ECM components and regulate immune cell
adhesion, migration, and polarization^[164]41. Recent studies have
indicated that COL6A3 can influence macrophage polarization toward an
M2-like phenotype, promoting an immunosuppressive
microenvironment^[165]52. Furthermore, the C-terminal fragment of
COL6A3, Endotrophin, has been implicated in modifying both tumor and
immune cell behavior, potentially facilitating macrophage recruitment
and functional shifts^[166]26.
In addition, COL6A3-mediated signaling has been explored as a potential
therapeutic target, with strategies such as blocking its interactions
or modulating its expression through antisense oligonucleotide therapy
currently being investigated in preclinical studies. These findings
further underscore the functional role of COL6A3 in shaping a
pro-tumorigenic microenvironment in recurrent tumors. In contrast,
primary tumors demonstrated enrichments of inflammatory macrophages and
monocytes, with the BMP5 signaling axis driving key cell-cell
interactions. Focusing on lymphocytes, we observed that the SPP1-CD44
and SPP1-ITGA4/ITGB1 signaling axes were the primary interactions
between immunogenic recurrent tumor cells and various lymphocytes,
which we speculate to inhibit T cell proliferation. Furthermore, HLA-E
crosstalk was particularly enriched between tumor cells with NKG2A + T
cells and NK cells, where NKG2A has been implicated in driving T cell
exhaustion in multiple tumor types.
From a therapeutic perspective, targeting COL6A3 and its downstream
effects offers several promising strategies. Blocking COL6A3-Integrin
α9β1 interactions could disrupt tumor-associated macrophage (TAM)
recruitment and polarization, potentially reducing immunosuppressive
signals within the TME. In addition, targeting COL6A3-mediated
β-catenin, TGF-β/Smad signaling pathways and endotrophin inhibition
also have been explored as a means to modulate fibrosis and immune
suppression, positioning it as a promising target in the context of
COL6A3-driven tumor progression. While direct inhibitors of COL6A3 have
not yet been developed, these findings suggest that targeting COL6A3,
either directly or through its interactions with integrins and
macrophages, could represent a viable approach for modulating the TME
and inhibiting tumor progression. Further research is warranted to
explore and validate these therapeutic avenues, which could ultimately
lead to the development of innovative treatments for COL6A3-driven
malignancies.
In summary, our study delineates the longitudinal evolutionary
trajectory of meningioma at single-cell resolution. Our results show
that recurrent meningiomas are characterized by distinct
transcriptional cellular states and tumor microenvironment compositions
compared to primary tumors. We also observed extensive diversity in
meningioma heterogeneity and evolution, with recurrent tumors evolving
toward a malignant and proliferative transcriptional trajectory.
Collectively, these findings advance our understanding of the
mechanisms underlying longitudinal meningioma biology and highlight
potential therapeutic targets for the treatment of high-risk recurrent
meningioma patients.
Methods
Sample collection and processing
With appropriate approval from the institutional review board (IRB),
meningioma specimens were collected from patients undergoing surgery at
Samsung Medical Center. The study protocol was sanctioned by our
institution’s ethical committees (Samsung Medical Center), and written
informed consent was obtained from all participants.
Single-nuclei RNA-sequencing data process and analysis
To process snRNA-seq, we utilized Cell Ranger (v.5.0.0) from 10x
Genomics with Count functionality for aligning reads to the prebuilt
GRCh38 genome reference v.2020-A (refdata-gex-GRCh38-2020-A). The
resulting gene-by-cell unique molecular identifier (UMI) count matrix
was then employed by the R package Seurat (v.4.0.5)^[167]53 for
subsequent analysis. Quality control measures were applied to remove
barcodes falling into specific categories: those likely representing
debris (expressing < 200 genes), potential doublets (expressing > 3700
genes), and cells showing signs of stress or apoptosis (with > 10%
mitochondrial gene expression). Following cell filtering, raw feature
counts were log-normalized, scaled, and underwent linear dimensional
reduction using principal component analysis (PCA). To correct for
batch effects across different samples, the Harmony (v.0.1.0) was
applied to the PCA embeddings^[168]54. Cell clusters were identified
based on a K-nearest neighbor (KNN) graph model computed on Euclidean
distances in PCA space, and visualization was performed using the
t-distributed stochastic neighbor embedding (tSNE) algorithm. Cell type
identities were assigned by performing differential gene expression
analysis within each cluster and annotating based on the expression of
characteristic markers.
Copy number alteration analysis
To detect large-scale chromosomal copy number variations (CNVs) using
snRNA-seq data, we employed inferCNV (v.1.6.0) with default parameters
recommended for 10x Genomics data^[169]55. This analysis was conducted
at the sample level using the raw counts matrix obtained after
post-quality-control filtering. Prior to running inferCNV across all
samples, we adjusted the bash environment with ‘ulimit -s unlimited’
and set specific options (expressions = 500000) within R.
Differential expression analysis
For identifying differentially expressed genes between distinct
clusters, we utilized the t.test function in Seurat, applying multiple
hypothesis correction via the Benjamini–Hochberg procedure. Genes with
an adjusted p-value < 0.05 and an absolute log2-fold change greater
than 1 were considered differentially expressed. Log2-fold changes were
computed based on the difference in log2-transformed mean counts
between groups.
RNA velocity analysis
To quantify spliced and unspliced unique molecular identifiers (UMIs)
per gene per cell, we employed the Python package velocyto^[170]56;
followed by subsequent analyses using scanpy and scvelo^[171]19. For
each cell, we computed moments (means and uncentered variances) of
normalized spliced and unspliced counts using scv.pp.moments. RNA
velocity estimation was performed with scv.tl.velocity using the
‘dynamical’ mode, and a velocity graph depicting transition
probabilities was constructed using scvelo.tl.velocity_graph.
Molecular subtype classification
Using the molecular gene set for each subgroup from the previous
study^[172]10, the ssGSEA method was applied to quantify the activity
of each molecular subtype within individual tumor cells. All genes were
ranked based on their expression levels, and enrichment scores for each
molecular subtype were computed by calculating a running sum statistic
that evaluates whether genes in each subtype gene set are enriched at
the top of a ranked list compared to genes not included. This enabled
quantification and comparison of the activity of the four molecular
subtypes across individual tumor cells. To mitigate potential technical
biases, molecular subtype scores were normalized to ensure that
differences in gene set size did not skew the results.
Two-dimensional representation of meningioma subtype states
Meningioma tumor cells were first separated into Immune/NF2 WT versus
Hypermetabolism/Proliferative by the sign of D, where D was calculated
as:
[MATH:
D=max(S
CImmune,SCNF2wt)−max(<
mrow>SCHy<
mi>per,S<
/mi>CPro
) :MATH]
1
SC[X] represents the single-sample gene set enrichment analysis
(ssGSEA)-derived score of the corresponding molecular subtype for a
given cell, measuring the enrichment level of the corresponding
biological signature within that cell. The D value was used to define
the y-axis of all cells.
Next, for Immune/NF2 WT cells (i.e., D > 0), the x-axis value was
defined as:
[MATH:
log2
(∣SC
Immune<
mo>−SCNF2wt∣+
1) :MATH]
2
For Hypermetabolism/Proliferative cells (i.e., D < 0), the x-axis was
defined as:
[MATH:
log2
(∣SC
Hyper−<
msub>SCPro∣+1)
:MATH]
3
The resulting differences have been used to plot cells based on the
corresponding quadrants representing the activity of each molecular
subtype.
Feature selection using correlation with gene expression and latent time
Using glmnet package^[173]57, elastic net regularization was
implemented with gene expression and latent time by alpha value = 0.5
and adopted a bootstrapping strategy 100 times to obtain the top 200
genes from regression analysis. Features were selected considering
their correlation values with gene expression, and latent time and the
number of appearances during the bootstrapping process.
Pathway analysis
Pathway enrichment analysis was performed using the clusterProfiler
package for functional analysis of differentially expressed genes in
recurrent compared to primary meningioma. Single-sample GSEA (ssGSEA)
was performed using the Easy single-cell analysis platform for
enrichment (escape, version 1.9.0)^[174]58. Enrichment scores were
calculated using the c2.all.v7.1.symbols.gmt curated gene set.
Survival Analysis Based on COL6A3 Expression Levels
To classify high and low expression levels of COL6A3, the maxstat R
package was applied to the bulk RNA sequencing dataset from an external
cohort of 110 meningioma patients^[175]31. The optimal cutoff was
determined using maxstat R package based on survival data. Survival
analysis was performed to compare overall survival (OS) between
patients with higher and lower COL6A3 expression using Kaplan-Meier
curves and log-rank tests.
Cell culture
The human meningioma cell lines IOMM-Lee and SF3061 were kindly
provided by Dr. Shin Hyuk Kang (Department of Neurosurgery, Korea
University Anam Hospital), who received them from the original
investigators^[176]21. IOMM-Lee was originally established by Dr. Wei
Hua Lee^[177]59, and SF3061 by Dr. Anitha Lal^[178]60. Cells were
regularly confirmed to be free of mycoplasma by PCR tests and
authenticated for matching short tandem repeat DNA profiles of the
original cell lines. The short tandem repeat profile of IOMM-Lee cell
line was 100% matched with the ATCC database, whereas SF3061, which is
not included in the ATCC database, was confirmed free of contamination
from other cell lines. Cells were cultured in Dulbecco’s Modified
Eagle’s Medium (DMEM; Welgene, South Korea) supplemented with 10% fetal
bovine serum (FBS; Avantor, PA, USA) and 1% penicillin-streptomycin
(Gibco, NY, USA). Cells were maintained at 37 °C in a humidified
incubator with 5% CO[2].
Reagents and transfection
Lipofectamine™ RNAiMAX Transfection Reagent (#13778), RNase A, TRIzol®
Reagent, and PowerUp SYBR™ Green Master Mix were obtained from
Invitrogen (CA, USA). siRNAs used included Silencer™ Select Negative
Control No.1 (#4390843) and Silencer™ Select COL6A3 siRNA (#S3318)
(Invitrogen). For cell viability assays, the Premix WST-1 Cell
Proliferation Assay System (TaKaRa, Japan) was used. Propidium iodide
(PI; Sigma-Aldrich, MO, USA) was dissolved in PBS and used for
staining. For cDNA synthesis, TOPscript™ RT DryMIX (Enzynomics, South
Korea) was utilized according to the manufacturer’s instructions.
siRNA Transfection
The cells were seeded (1.5 × 10^3 or 3 × 10^4 cells per well,
respectively) in a 96-well or 12-well flat-bottom plate and incubated
overnight. They were then transfected with siRNA using Lipofectamine
RNAiMAX Transfection Reagent.
Cell proliferation assay
Seventy-two hours after transfection, 10 μL of Premix WST-1 (TaKaRa)
was added to each well and the plate was incubated for 2 h. The
absorbance was measured at 420 nm using a VersaMax Microplate reader
(Molecular Devices)
Cell cycle assay
Twenty-four hours after transfection, the cells were cultured in
serum-free medium for starvation for an additional 24 h. Subsequently,
the medium was replaced with serum-containing medium for stimulation.
After 6 h, the cells were harvested and fixed with ethanol at 4 °C
overnight. Following fixation, the cells were washed with PBS and
treated with RNase A (100 μg/mL) at 37 °C for 15 min. The cells were
then stained with PI (50 μg/mL). Finally, PI-stained cells were
analyzed using a BD LSRFortessa™ X-20 (BD Biosciences) and FlowJo^TM
software version 10.10.0.
Quantitative real-time PCR (qRT-PCR)
Forty-eight hours after transfection, the cells were harvested using
TRIzol® Reagent. Total RNA was isolated following the manufacturer’s
user guide (Invitrogen). cDNA was synthesized using TOPscript™ RT
DryMIX (enzynomics). Quantitative real-time PCR (qRT-PCR) assays were
performed on a QuantStudio^TM 3 system using SYBR Green master mix. The
primer sequences for GAPDH were as follows: sense 5’-
CATGTTCGTCATGGGTGAACCA-3’ and antisense 5’-
AGTGATGGCATGGACTGTGGTCAT-3’. The primer sequences for COL6A3 were as
follows: sense 5’- TCTCTTAAAATCAGTGCACAACG-3’ and antisense 5’-
AACTCTTTCAACAGAGGGAAGC-3’. The results were analyzed using the 2^-ΔΔCt
method.
Histology and immunohistochemistry
Formalin-fixed, paraffin-embedded tumor sections (4 μm) were stained
with hematoxylin and eosin (H&E) and processed for immunohistochemistry
(IHC). Antigen retrieval was performed using either Tris-EDTA buffer
(60 °C, overnight) or 0.1% trypsin (37 °C, 30 min). Tissue sections
were blocked with 3% hydrogen peroxide and 5% bovine serum albumin
(BSA) in PBS. The primary antibodies used were rabbit anti-Collagen VI
alpha 3 (polyclonal, MYBioSource; catalog number MBS9612481; dilution
1:100), anti-CD163 (clone EPR14643-36, Abcam; ab189915; RRID:
AB_3493124; dilution 1:1000), and anti-C1QA (clone [179]EPR14634,
Abcam; ab189922; RRID: AB_2894866; dilution 1:1000). A biotinylated
anti-rabbit secondary antibody (Jackson ImmunoResearch; dilution
1:1000) was applied, followed by Streptavidin-HRP (Invitrogen; dilution
1:1000) and visualization with 3,3′-diaminobenzidine (DAB; Abcam).
Finally, sections were counterstained with hematoxylin. Representative
images were acquired using a Zeiss Axioscan Z1 Slide Scanner, and
regions of interest (ROIs) of consistent size were selected and cropped
using ZEN software for downstream analysis.
Cell‒cell communication analysis
To analyze cell-to-cell communication involving tumor cells and other
cell types, we utilized the R package CellChat (v.1.6.1)^[180]61 with
the ligand‒receptor interaction database “CellChatDB.human”. The
default preprocessing steps were applied, and functions
computeCommunProb and computeCommunProbPathway were used to infer
networks for individual ligand‒receptor pairs and signaling pathways.
Cell type proportion analysis using miloR
Differences in cell proportions between primary and recurrent samples
were assessed using the miloR (v.1.9.1) package^[181]62 following a
standard workflow. This involved building graphs and defining
neighborhoods (similar cells identified by k-nearest neighbors (kNN))
using the buildGraph and makeNhoods functions with parameters
prop = 0.1, k = 100, d = 30. Parameters were fine-tuned based on manual
inspection of neighborhood sizes to optimize informative distributions
for downstream analyses. Differential abundance within each
neighborhood was evaluated using the testNhoods function with false
discovery rate (FDR) adjustment based on neighborhood distances
calculated by calcNhoodDistance.
Statistics & reproducibility
Statistical analyses were conducted using R software (version 4.0.5).
Differences in categorical variables (e.g., cell type proportions) were
assessed using Fisher’s exact test, while continuous variables were
compared using Student’s t-test or Wilcoxon rank-sum test, as
appropriate. All hypothesis tests were two-sided, and statistical
significance was set at p < 0.05. No statistical method was used to
predetermine sample size, and the experiments were not randomized.
Reporting summary
Further information on research design is available in the [182]Nature
Portfolio Reporting Summary linked to this article.
Supplementary information
[183]Supplementary Information^ (58.5MB, pdf)
[184]Reporting Summary^ (89.1KB, pdf)
[185]Transparent Peer Review file^ (5.5MB, pdf)
Source data
[186]Source Data^ (23.3MB, xlsx)
Acknowledgements