Abstract
Background
Intratumoral heterogeneity (ITH) and tumor microenvironment (TME) of
pancreatic ductal adenocarcinoma (PDAC) play important roles in tumor
evolution and patient outcomes. However, the precise characterization
of diverse cell populations and their crosstalk associated with PDAC
progression and metastasis is still challenging.
Methods
We performed single-cell RNA sequencing (scRNA-seq) of treatment-naïve
primary PDAC samples with and without paired liver metastasis samples
to understand the interplay between ITH and TME in the PDAC evolution
and its clinical associations.
Results
scRNA-seq analysis revealed that even a small proportion (22%) of
basal-like malignant ductal cells could lead to poor chemotherapy
response and patient survival and that epithelial-mesenchymal
transition programs were largely subtype-specific. The clonal
homogeneity significantly increased with more prevalent and pronounced
copy number gains of oncogenes, such as KRAS and ETV1, and losses of
tumor suppressor genes, such as SMAD2 and MAP2K4, along PDAC
progression and metastasis. Moreover, diverse immune cell populations,
including naïve SELL^hi regulatory T cells (Tregs) and activated
TIGIT^hi Tregs, contributed to shaping immunosuppressive TMEs of PDAC
through cellular interactions with malignant ductal cells in PDAC
evolution. Importantly, the proportion of basal-like ductal cells
negatively correlated with that of immunoreactive cell populations,
such as cytotoxic T cells, but positively correlated with that of
immunosuppressive cell populations, such as Tregs.
Conclusion
We uncover that the proportion of basal-like subtype is a
key determinant for chemotherapy response and patient outcome, and that
PDAC clonally evolves with subtype-specific dosage changes of
cancer-associated genes by forming immunosuppressive microenvironments
in its progression and metastasis.
Supplementary Information
The online version contains supplementary material available at
10.1186/s12943-024-02003-0.
Keywords: Pancreatic ductal adenocarcinoma, Liver metastasis,
Single-cell RNA-sequencing, Intratumoral heterogeneity, Tumor
microenvironment
Introduction
Pancreatic ductal adenocarcinoma (PDAC) is one of the most lethal
diseases [[55]1]. Most PDAC patients present with nonspecific symptoms
at an advanced stage that is not amenable to curative surgery.
Moreover, PDAC has clinical characteristics of early metastasis and
resistance to therapy [[56]2]. However, molecular mechanisms underlying
PDAC progression, metastasis, and therapy response are still poorly
understood.
With recent advances in single-cell RNA sequencing (scRNA-seq)
technologies, it is now possible to characterize comprehensively the
cellular composition and the transcriptomic landscape of the tumor
microenvironment (TME) at the single-cell level. Therefore, scRNA-seq
can be useful in studying the complex nature of intratumoral
heterogeneity (ITH) and TME in PDAC for establishing novel therapeutic
strategies. A few scRNA-seq studies of PDAC have been published, but
the majority of these works are limited to either primary tumors
[[57]3–[58]5] or liver metastases [[59]6]. Furthermore, surgical
specimens are usually acquired after portal vein dissection and
vascular clamping, which would cause ischemic damage and autolysis of
pancreatic cells. However, endoscopic ultrasound-guided fine needle
biopsy (EUS-FNB) enables acquiring relatively undamaged and fresh
tissue under the original physiologic status of the patient [[60]7].
Recently, Zhang et al. reported that CEACAM5^+/CEACAM6^+ ductal cells
are associated with poor prognosis by analyzing three cases of primary
PDAC and matched liver metastases. However, these datasets offer
restricted insights into PDAC evolution due to their focus on only the
advanced tumor stage and a relatively small sample size [[61]8].
Here, we performed scRNA-seq analysis of treatment-naïve 21 primary
PDAC samples obtained by EUS-FNB and 7 matched liver metastasis samples
acquired by percutaneous biopsy at the time of diagnosis to better
understand the interplay between ITH and TME in tumor evolution and its
clinical relevance in the treatment of PDAC.
Results
Patient characteristics
This study enrolled 21 treatment-naïve patients (Table S[62]1). The
median age was 61 (50–73 years), and 13 patients (62%) were females.
Tumor clinical stages (the 8th AJCC) were 6 (29%) at stage III and 15
(71%) at stage IV. Among the 15 patients with stage IV disease, 13 had
metastasis to the liver, and two had metastasis to the bone or the
lymph node but not the liver. The median overall survival (OS) was
9.7 months, ranging from 0.6 to 47.8 months.
Single-cell transcriptional landscape of primary PDACs and matched liver
metastases
From these 21 patients, we obtained the following samples; 1) primary
PDAC without metastasis (Pm0, N = 6), 2) primary PDAC with metastasis
(Pm1, N = 15), 3) liver metastasis matched with primary PDAC (Lm,
N = 7), and 4) adjacent normal pancreas (Pn, N = 5) (Fig. [63]1A, Table
S[64]2). We divided cells into 26 clusters and identified seven major
cell types (Fig. S[65]1A and B) ([66]Supplementary Methods). We did not
observe a batch effect, although the ductal cell clusters showed
evident patient-, rather than origin-, specific gene expression
profiles, which is typical in tumor cells with patient-specific copy
number variations (CNVs) [[67]5, [68]6] (Fig. S[69]1C-E). When compared
to previous scRNA-seq studies of PDAC, T cells were relatively enriched
while fibroblasts were somewhat depleted in our data. This may be due
to the way of acquiring tissues with EUS-FNB in our study, while most
of the previous works used surgical resection.
Fig. 1.
[70]Fig. 1
[71]Open in a new tab
ScRNA-seq analysis of PDAC subtypes and their clinical relevance. A
Schematic of the experimental design. ScRNA-seq was performed on PDAC
samples from 21 patients, including non-metastatic PDACs (N = 6),
metastatic PDACs (N = 15), and matched liver metastases (N = 7). B
Heatmap showing the expression of signature genes for NMF subtypes.
Each column in the heatmap corresponds to one cell and each row of the
heatmap corresponds to a signature gene of four NMF subtypes. Origin,
patient, NMF subtype, and previously reported PDAC classification
schemes for each cell are shown at the top of the plot, and the results
of GSEA for each signature gene set are shown on the right side of the
plot. C Pathway enrichment analysis focusing on origin-specific
differences within classical and basal-like subtypes. D and E
Kaplan–Meier overall survival curves for PDAC patients based on the
fraction of basal-like subtype in the deconvoluted TCGA PDAC RNA-seq
dataset (D), and in their primary PDAC in our dataset (SMC cohort) (E).
F Forest plot showing the estimated hazard ratios for the
clinicopathologic parameters and the proportions of NMF subtypes by
multivariate Cox regression analysis of combined scRNA-seq data from
our cohort and the two previously published PDAC cohorts. Data are
presented as hazard ratio ± 95% confidence interval. G Waterfall plot
showing the best percentage change in the sum of the target lesions
according to the RECIST v1.1. Each bar indicates a study sample, and
the sample is divided into two groups: those with basal-like proportion
above 22% (red) and those below (cyan). H and I The proportion of PDAC
NMF subtypes and CT scan images before and after chemotherapy of PDAC
patients PB2341 (H) and PB2311 (I). J Boxplot showing the distribution
of mean CNV correlation coefficients among malignant ductal cells
within origins (two-sided Wilcoxon rank sum test: *P < 0.05,
**P < 0.01, ***P < 0.001). K and L Hierarchical clustering of CNV
profiles in individual patients PB2155 (K) and PB2191 (L). M and
N Unsupervised transcriptional trajectories of ductal cells in
individual patients PB2155 (M) and PB2191 (N) colored by sample origin.
Trajectory directions were indicated by arrows. O and P Dots on
trajectory projections (left) were colored by copy number scores at the
cellular level and overlaid with contour plots of cells with the
strongest copy number variation for known cancer-associated genes in
individual patients PB2155 (O) and PB2191 (P). Violin plots (right)
showed copy number scores of genes by origin (two-sided Wilcoxon rank
sum test: *P < 0.05, **P < 0.01, ***P < 0.001)
Stratification of PDAC subtypes
Sub-clustering analysis classified ductal cells into 21 subclusters,
each of which was largely patient-specific (Fig. S[72]2A). To identify
distinct transcriptional programs, we applied consensus non-negative
matrix factorization (cNMF) on the gene expression profiles of ductal
cells. Among the four NMF subtypes determined by stability and error
(Fig. S[73]2B), three NMF subtypes were well matched with the
previously reported PDAC ([74]Supplementary Results). Therefore, we
hereafter designated NMF-1 as ‘cycling’, NMF-2 as ‘classical’, NMF-3 as
‘normal-like’, and NMF-4 as ‘basal-like’.
Since the ‘epithelial–mesenchymal transition (EMT)’ pathway was
involved in both classical and basal-like subtypes (Fig. [75]1B), we
further investigated different EMT mechanisms between the two subtypes
through pseudobulk-based differential gene expression analysis and gene
set enrichment analysis (GSEA) among origins. Interestingly, the 'EMT'
pathway was significantly activated in Pm1 of both classical and
basal-like subtypes (Fig. [76]1C); however, actual genes involved in
the 'EMT' pathway were different between the two subtypes, implying
subtype-specific EMT programs might exist in PDAC (Fig. S[77]3C). To
interrogate factors potentially contributing to the subtype-specific
EMT programs, we analyzed underlying transcription factors (TFs)
associated with the distinct subtype-specific EMT genes in Pm1. As
shown in the Fig. S[78]3D and E, 38 and 62 TFs were identified as
Pm1-specifically activated regulons in the classical and basal-like
PDAC subtypes, respectively. Among the Pm1-specifically activated TFs
in the classical subtype, we found that TFDP1 and CUX1 were known to
respectively regulate IGFBP2 and MFAP5 that were detected as the
classical-specific EMT-associated differentially expressed genes (DEGs)
in Pm1 (Fig. S[79]3D). We also discovered that 16 out 62
Pm1-specifically activated TFs in the basal-like subtype are involved
in regulating all the basal-like-specific EMT-associated DEGs in Pm1
except COL1A2 (Fig. S[80]3E). Interestingly, the two Pm1-specifically
activated TFs in the classical subtype and the 16 Pm1-specifically
activated TFs in the basal-like subtype did not overlap with each
other.
PDAC subtypes and clinical outcomes
We identified that malignant ductal cells with different subtypes
coexist in individual samples and that the proportion of each subtype
varied across the samples even within the same patient (Fig. S[81]3F).
Interestingly, among seven patients having both matched primary PDAC
and liver metastasis samples, five patients (PB2191, PB2264, PB2349,
PB2409, and PB2410) shared common predominant subtypes, but two
patients (PB2155 and PB2311) exhibited different predominant subtypes
(Fig. S[82]3F).
As the basal-like subtype of PDAC is known to be associated with a poor
prognosis [[83]9], we interrogated which proportion of basal-like
subtype is associated with patients’ survival. When we tested a range
of basal-like cell proportions (10% ~ 35%) while maintaining the number
of samples in each group greater than 10% of total samples based on the
external PDAC scRNA-seq datasets (WashU [[84]4] and MGH [[85]10],
N = 25), the lowest cutoff for basal-like cell proportion showing a
statistical association with survival was 22% (P = 0.024, Fig.
S[86]4A). In addition, we deconvoluted bulk RNA-seq data from TCGA PDAC
cohort [[87]11] (N = 148) using our scRNA-seq data and estimated
cellular fraction of the four PDAC subtypes in each TCGA PDAC sample.
When we scanned the proportions of basal-like cells and their
statistical association with survival based on the deconvoluted TCGA
PDAC data, 22% of basal-like cell fraction was the lowest threshold
showing a statistical association with survival (P = 0.033, Fig. [88]1D
and Fig. S[89]4A). This result was also consistent in our cohort (SMC,
N = 18) (P = 0.023, Fig. [90]1E). To estimate the prognostic relevance
of NMF subtypes, we also performed multivariate Cox regression analysis
for OS with age, sex, grade, and proportion of NMF subtypes in our data
combined with the two previously published scRNA-seq data [[91]4,
[92]10]. Only the proportion of basal-like was significantly associated
with shorter OS (Hazard ratio, 24.9; 95% CI, 2.02–310; P = 0.012;
Fig. [93]1F).
Evaluation of treatment response according to PDAC subtypes
We then explored the association of PDAC subtypes with treatment
responses (Fig. [94]1G-I, Fig. S[95]4D-G). In line with the result
above, the proportion of basal-like was a key determinant for
chemotherapy response. The mean best percentage change in the sum of
the target lesions according to the RECIST v1.1 increased by 39% in the
group contains more than 22% of basal-like subtype, and decreased by
34% in the group with less than or equal to 22% of basal-like subtype
(Fig. [96]1G). The change in the sum of the target lesions between the
two groups was statistically significant (P = 0.0002, Fig. S[97]4B). In
addition, the proportions of basal-like showed positive correlations
with changes in tumor dimension (r = 0.73, P = 4.9
[MATH: × :MATH]
10^–5) (Fig. S[98]4C).
Figure [99]1H, I and Fig. S[100]4D to G showed the 1^st line
chemotherapeutic response according to the proportion of PDAC subtypes.
For example, patient PB2341 having primary PDAC with mixed subtypes of
classical (56%) and basal-like (36%) did not respond to the first four
cycles of gemcitabine plus nab-paclitaxel (GnP) followed by aggressive
progression, and died with an OS of 5.3 months (Fig. [101]1H). In
addition, PB2256 with a high proportion (79%) of basal-like in primary
PDAC showed rapid progression of primary mass after four cycles of
FOLFIRINOX and poor response to the subsequent GnP treatment as well,
with progressive disease (PD) after three cycles (Fig. S[102]4D).
Interestingly, in the case of PB2311 with different subtype
compositions between primary PDAC (49% of basal-like and 42% of
classical) and its liver metastasis (81% of classical), the primary
tumor mass increased whereas the liver metastasis mass decreased after
FOLFIRINOX treatment (Fig. [103]1I). On the other hand, patient PB2366
with normal-like predominant known to be excellent prognosis (93%)
primary PDAC had partial response (PR) after four cycles of FOLFIRINOX
treatment as the first-line chemotherapy, and the tumor was down-staged
to resectable status. The patient received pylorus-preserving
pancreatoduodenectomy with adjuvant FOLFIRINOX chemotherapy and has now
been followed up for 44.9 months with no evidence of disease (Fig.
S[104]4E). Patient PB2032 with mixed subtypes of normal-like (58%) and
classical (41%) in primary PDAC showed the best response as PR with
palliative GnP treatment and was still alive with 45.6 months of OS
(Fig. S[105]4F). Furthermore, patient PB2191 with more than 90% of
classical in both primary PDAC and liver metastasis had favorable
responses to chemotherapeutic drugs in both sites with 17.8 months of
OS (Fig. S[106]4G).
Dynamics of clonal evolution during PDAC progression
To understand clonal heterogeneity and evolution in the progression of
primary PDAC to liver metastases, we performed CNV analysis of ductal
cells from both primary PDAC and liver metastases. The CNV profiles
were highly patient-specific but still largely concordant with those of
TCGA PDAC cohort (Fig. S[107]5A). We analyzed the average CNV
correlation coefficients among malignant ductal cells in each sample to
measure the level of clonal heterogeneity in the tumor (Fig. [108]1J).
Interestingly, the average CNV correlation coefficients were lowest in
Pm0-derived malignant ductal cells and highest in Lm-derived ones, and
the differences were significant. The level of clonal heterogeneity of
Pm1-derived malignant ductal cells was in between Pm0 and Lm
(Fig. [109]1J, Fig. S[110]5B). This result suggests that the clonal
heterogeneity decreases as the tumor progresses and metastasizes to the
liver.
To understand the clonal evolution of primary PDAC to liver metastasis
in individual patients, we performed a hierarchical clustering analysis
of ductal cells based on their CNV profiles for each patient
(Fig. [111]1K and L, Fig. S[112]5C and D). Overall, CNV profiles of
primary PDACs and their matched liver metastases were generally
concordant to each other. However, several primary PDAC- or liver
metastasis-dominant CNV events were identified. We further investigated
copy number changes along the tumor evolution by performing trajectory
analysis of ductal cells in each patient. The lineage differentiations
from Pn-derived ductal cells into Pm1- and eventually Lm-derived
malignant ductal cells were well reconstructed (Fig. [113]1M and N,
Fig. S[114]5E and F). Notably, copy number gains of oncogenes or losses
of tumor suppressor genes gradually become prevalent and pronounced
along with PDAC evolution. In the cases where basal-like was
predominant in liver metastases, the CNV score of KRAS showed a
significantly positive correlation with pseudotime (Fig. [115]1O,
PB2155: r = 0.26, P = 1.9
[MATH: × :MATH]
10^–6; Fig. S[116]5G, PB2349: r = 0.20, P = 3.0
[MATH: × :MATH]
10^–15). Furthermore, malignant ductal cells with the top 10% CNV score
of KRAS were mostly Lm-derived cells, and the CNV score of KRAS was
also significantly higher in Lm than in Pm1 (Fig. [117]1O, Fig.
S[118]5G). In contrast, CNV score of tumor suppressor genes such as
SMAD2 and MAP2K4 showed a negative correlation with pseudotime, and the
copy number losses of SMAD2 and MAP2K4 were more prominent in the Lm of
PB2155 and PB2349, respectively (Fig. [119]1O, Fig. S[120]5G). In the
classical dominant cases, the CNV score of ETV1 had a significant
positive correlation with pseudotime (Fig. [121]1P, PB2191: r = 0.20,
P = 3.9
[MATH: × :MATH]
10^–7; Fig. S[122]5H, PB2264: r = 0.11, P = 4.1
[MATH: × :MATH]
10^–8), but the CNV score of KRAS did not significantly correlate with
pseudotime. The malignant ductal cells with the top 10% CNV score of
ETV1 were mostly observed in Lm, and the CNV score of ETV1 was
significantly higher in Lm than in Pm1 (Fig. [123]1P, Fig. S[124]5H).
The CNV scores of two other oncogenes NFE2L2 and PIK3CB also showed
positive correlations with pseudotime, and their CNV scores were
significantly higher in Lm than in Pm1 of PB2191 and PB2264,
respectively (Fig. [125]1P, Fig. S[126]5H).
To further interrogate origin- and subtype-specific CNVs, single-cell
CNV profiles were merged into the sample level and the frequencies of
CNV events were measured across samples. Copy number gains of
chromosome 2q31, 2q32, and 8q24 and losses of 6q21, 6q22, 18q12, 18q21,
and 18q22 occurred more frequently in Pm1 compared to Pm0. Due to these
differential copy number alterations between Pm1 and Pm0,
metastasis-associated genes such as NFE2L2 and EXT1 were amplified
while metastasis suppressor genes such as FOXO3, GOPC, PTPRK, SETBP1,
and SMAD2 were deleted in Pm1 (Fig. S[127]6A). Interestingly, Lm showed
significantly more frequent copy number gain of the chromosomal region
containing PPFIBP1 (12p11) that was known to be associated with tumor
development, progression, and metastasis of PDAC than Pm1 (Fig.
S[128]6B). In addition, copy number gain of chromosome 12q11-12
containing KRAS occurred more frequently in the basal-like than in the
classical (Fig. S[129]6C), which was consistent with previous reports
that KRAS amplification was more prominent in the basal-like subtype
[[130]6, [131]12]. However, copy number gain of chromosome 7p21
spanning ETV1 that was reported to promote pancreatic cancer metastasis
was more frequently observed in the classical than in the basal-like.
Niche and subtype-specific characteristics of T and NK cells
We sub-clustered T/NK cells to analyze their functional characteristics
in primary PDACs and liver metastases and identified a total of 24
subclusters (Fig. S[132]7A), which were subsequently classified by the
expression profiles of canonical marker genes (Fig. S[133]7B).
Differential gene expression analysis was also conducted to further
characterize each subcluster (Fig. S[134]7C and D).
When comparing Pm1 with Pm0, a naïve/resting regulatory T cell (Treg)
subcluster, Treg-SELL, and an unstimulated natural killer (NK) cell
subcluster, NK-XCL2, were notably enriched in Pm1. We also found that
the proportions of helper T cell (Th)-GRP183 and NK-KLRC2 with
antitumor properties were significantly reduced in Lm compared to in
Pm0 and Pm1 and to in Pm1, respectively. In contrast, exhausted T
cell (Tex)-LAG3 and Treg-TIGIT with immunosuppressive characteristics
were remarkably enriched in Lm compared to in Pm1 and to in Pm0 and
Pm1, respectively. A dysfunctional NK cell subcluster, NK-KLRC1, was
also marginally enriched in Lm (Fig. [135]2A, Fig. S[136]8A). These
patterns were also identified in individual patients with primary PDAC
and matched liver metastases (Fig. [137]2B, Fig. S[138]8B and C). The
trajectory analysis of CD4^+ T cells, CD8^+ T cells, and NK cells
confirmed that their regulatory natures and dysfunctional
characteristics gradually became prominent along PDAC progression and
metastasis (Figs. S[139]9 and S[140]10).
Fig. 2.
[141]Fig. 2
[142]Open in a new tab
The interplay between ITH and TME in the primary PDACs and matched
liver metastases. A Box plots indicating the percentage differences in
T cell subclusters among origins (two-sided Wilcoxon rank sum test:
*P < 0.05, **P < 0.01, ***P < 0.001). B Area plots displaying the
changes in T cell subcluster composition by origin for each patient. C
Dot plots illustrating ligand-receptor interactions between malignant
ductal cells and Tregs. The size of a circle indicates an interaction
score, and the color of a circle represents the origin. D Multiplex
immunohistochemistry (IHC) showing the interaction between ICAM1
(magenta)- or IGF2R (orange)-expressing FOXP3^+ (green) Tregs (arrows)
and AREG (cyan)- or IGF2 (yellow)-expressing cytokeratin (CK)^+ (red)
tumor cells. Nuclei are counterstained with DAPI (blue). E and F
Scatter plot displaying the correlation between the fraction of
basal-like in ductal cells and the fraction of cytotoxic T cells in T
cells (E), and between the expression level of S100A9 in ductal cells
and the fraction of cytotoxic T cells in T cells (F). G and H Pearson
correlation between the proportion of basal-like among ductal cells and
the proportion of Tregs among the T cell population (G), and between
the expression level of S100A9 in ductal cells and the fraction of
Tregs among T cells (H). I and J Mapping of major cell types (I) and
NMF subtypes (J) to spatial transcription spots from treatment naïve
PDAC patient published by Zhou et al. using a robust cell type
decomposition (RCTD) method. K The spots on the spatial transcriptome
slide were colored by NMF subtypes and overlaid with contour plots of
Treg enriched spots. L and M Multiplex IHC showing the expression of
S100A9 (yellow) and the distribution of T cells in basal-like dominant
(L) and classical dominant (M) PDAC tissues. CD8 (green) for cytotoxic
T cells, FOXP3 (red) for Tregs, CK (white) for ductal cancer cells,
S100A2 (magenta) for basal-like ductal cells and DAPI (blue) for nuclei
were co-stained. Scale bar, 50 μm. N Box plots indicating the
percentage differences in myeloid subclusters among origins (two-sided
Wilcoxon rank sum test: *P < 0.05, **P < 0.01, ***P < 0.001). Samples
from the same patients were connected by solid lines. O Area plots
showing the change in the composition of the myeloid subclusters by
origin for each patient. P and Q Scatter plot displaying the Pearson
correlation between the fraction of basal-like in ductal cells and the
fraction of Mono-FCN1 (P) and Mp-TGFBI (Q)
According to the above results, we hypothesized that the incremental
activation of Tregs in PDAC evolution was potentially attributed to
their interaction with malignant ductal cells, and investigated
cellular interactions between ductal cells and Tregs. We discovered
that Treg-activating intercellular interactions gradually established
as PDAC progressed and metastasized into the liver. For example,
Treg-stabilizing LGALS9-CD44 interactions between ductal cells and
Tregs were identified in all three origins (Fig. [143]2C). However,
FOXP3-inducing TNFSF9-TNFRSF9 interactions were observed in Pm1 and Lm,
but not in Pm0. Interestingly, an IGF2-IGF2R interaction known to
promote Treg proliferation was only observed in Lm, and another
Treg-enhancing AREG-ICAM1 interaction between ductal cells and
Treg-TIGIT was also identified uniquely from Lm (Fig. [144]2C). These
two liver metastases-specific immunosuppressive interactions between
ductal cells and Tregs were further validated using multiplex
immunohistochemistry (mIHC) (Fig. [145]2D, Fig. S[146]11-[147]13).
To interrogate the association between PDAC subtypes and the immune
environment, we analyzed correlations between the proportions of ductal
and T cell subtypes. The proportion of basal-like was negatively
correlated with the proportion of cytotoxic T cells (Fig. [148]2E), as
previously reported by Raghavan et al. [[149]6] and Hwang et al.
[[150]4]. Especially, S100A9, a basal-like signature gene known to
suppress T cell proliferation by activating myeloid-derived suppressor
cells, showed significant negative correlation (Pearson's correlation
r = -0.63; P = 0.001) between its expression level and the proportion
of cytotoxic T cells (Fig. [151]2F). In contrast, the proportion of
Tregs was positively correlated with the proportion of basal-like
(Pearson's correlation r = 0.42; P = 0.048, Fig. [152]2G), and the
expression level of S100A9 showed a positive relationship with the
fraction of Tregs in T cells (Pearson's correlation r = 0.43;
P = 0.039, Fig. [153]2H). To validate these results, we performed a
deconvolution analysis of PDAC spatial transcriptomic data published by
Zhou et al. [[154]10] using our scRNA-data as a reference
(Fig. [155]2I-K). Tregs were largely enriched in the regions where
basal-like was predominant, and a significant positive correlation
between the proportions of basal-like and Tregs in spatial
transcriptomic spots was identified (Pearson’s correlation r = 0.135;
P < 0.001). The distribution of CD8^+ T cells and Tregs around
basal-like ductal cells expressing S100A2 and S100A9 were further
confirmed by mIHC (Fig. [156]2L and M, Fig. S[157]14).
TGFBI^hi macrophage shapes an immunosuppressive environment in liver
metastasis
Fourteen distinct clusters of myeloid cells were identified (Fig.
S[158]15A), and each cell type were specified with previously described
markers for myeloid (Fig. S[159]15B). Pro-inflammatory monocyte
(mono)-IL1B and anti-inflammatory mono-CDKN1C subclusters were
significantly increased in Pm1 while Mono-CCR2, macrophage (Mp)-C1QB,
Mp-LY6E, and Mp-TGFBI subclusters involved in tumorigenesis and
immunosuppression were enriched in Lm compared to Pm1 (Fig. [160]2N and
O, Fig. S[161]15C and D).
Furthermore, the proportion of mono-FCN1, which is a classical monocyte
involved in the initial inflammatory response, was inversely correlated
with the proportion of basal-like malignant ductal cells
(Fig. [162]2P). In contrast, the fraction of Mp-TGFBI with
immunosuppressive properties was positively correlated with the
fraction of basal-like cells (Fig. [163]2Q). Among 59 basal-like
signature genes in ductal cells, 32 genes were negatively associated
with Mono-FCN1 and 32 genes such as immunosuppression-related HCAR2 and
CTHRC1 were positively correlated with Mp-TGFBI (Fig. S[164]15E).
Discussion
We have explored the landscape of diverse cellular populations in
primary PDACs and their matched liver metastases to comprehensively
understand the molecular mechanisms of tumor progression, metastasis,
and treatment response in PDAC. Our work demonstrates the complex
nature of ITH and TME of PDAC at the single-cell level and its
association with chemotherapy responses and clinical outcomes. We also
discovered that TME becomes more immunosuppressive and that the clonal
heterogeneity of tumor decreases as PDAC progresses and metastasizes
into the liver.
Several studies have investigated the characteristics of PDAC subtypes
and their association with patient prognosis using bulk RNA-sequencing
data [[165]9, [166]12]. However, this approach has a fundamental
limitation in that only dominant subtypes can be identified, making it
difficult to interrogate the association between heterogeneous cellular
compositions of PDAC subtypes and clinical prognosis. In this study, we
revealed that even a small proportion (~ 22%) of basal-like malignant
ductal cells has a detrimental effect on patients' survival and the
overall response rate to chemotherapy, highlighting the advantages of
scRNA-seq in estimating the exact cellular proportion of specific
subtypes, which is more informative than the prevailing subtypes for
predicting a patient's prognosis.
EMT has been reported to be predominantly activated in the basal-like
subtype in PDAC [[167]12]. However, when comparing primary PDACs and
liver metastases separately by PDAC subtypes in our scRNA-seq data, the
EMT pathway was notably enriched in both classical and basal-like
ductal cells of metastatic primary PDACs, although the EMT-related
genes were largely subtype-specific. Interestingly, we discovered that
TFs known to regulate the subtype-specific EMT genes were
also activated in a subtype-specific fashion, which implies that gene
regulatory networks controlling EMT programs may be differently wired
in the two PDAC subtypes, resulting in the activation of
the subtype-specific EMT genes and indicating the need to develop
personalized treatment strategies considering the ITH of PDAC.
Furthermore, we revealed that clonal homogeneity of tumor cells
increases with PDAC progression and metastasis from the CNV analysis.
Although primary PDACs and liver metastases shared most of CNVs, CNVs
unique to either primary PDACs or liver metastases were also observed
at a single-cell level. This result confirms that metastatic tumor
cells evolve from one of clonal groups within the primary PDAC and that
continuous clonal expansion occurs in both the primary and metastatic
tissues even after metastasis [[168]13]. Our analysis further
identified CNV events conferring metastatic advantages. Notably, the
copy numbers of cancer-associated genes such as KRAS, ETV1, and SMAD2,
displayed a subtype-specific increase in prevalence and magnitude along
PDAC progression, highlighting the importance of matched primary and
metastatic samples for uncovering such patterns. In particular, our
findings are in line with the assertion by Mueller et al. that KRAS
dosage is an important contributor to PDAC progression and metastasis
[[169]14]. As previously reported, KRAS amplification was more
prominent in the basal-like subtype than the classical subtype in our
data [[170]12]. However, in the classical subtype, we found that ETV1
dosage exhibits a trend similar to KRAS dosage in the basal-like
subtype, with ETV1 gain being more widespread and pronounced along PDAC
progression and metastasis. These findings provide new avenues for
understanding the complex mechanisms of PDAC progression and metastasis
in a subtype-specific manner.
Immunosuppressive Tregs are accumulated in PDAC and pancreatic
intraepithelial neoplasia and related to metastasis and poor prognosis
[[171]3]. However, the impact of Tregs in PDAC progression and
metastasis still remains elusive. We identified two Treg subclusters,
Treg-SELL and Treg-TIGIT, significantly enriched in primary PDACs and
liver metastases, respectively. Treg-TIGIT is terminally differentiated
or activated Treg while Treg-SELL is in a naïve or resting state
[[172]15]. Notably, Treg-activating cellular interactions
were established stepwise during PDAC progression and metastasis. These
findings indicate that Tregs play an essential role in shaping the
immunosuppressive TME of both primary PDAC and liver metastasis and
that their functional characteristics become more regulatory as PDAC
progresses.
CD8^+ and CD4^+ T cells are depleted while C1QC^+ macrophages are
predominant in the TME of the basal-like subtype-dominant cases
[[173]6]. We confirmed that the proportion of CD8^+ T cells is
negatively correlated with that of basal-like ductal cells. However,
the proportions of CD4^+ T cells and C1QC^+ macrophages did not
correlate significantly with the proportions of basal-like subtype in
our scRNA-seq data. Instead, the proportions of C1QC^+ macrophages were
significantly more enriched in liver metastases than in primary PDACs.
Moreover, the proportions of classical monocytes (Mono-FCN1) were
negatively correlated with the proportions of basal-like ductal
subtype. In contrast, the proportions of Tregs and immunosuppressive
macrophages (Mp-TGFBI) were positively correlated with the basal-like
ductal subtype. Therefore, these suggest that the proportions of
basal-like subtype tend to be associated with more immunosuppressive
environments compared to the other subtypes.
Conclusions
Here, we provide a comprehensive catalog of ITH and TME from
non-metastatic PDAC and metastatic PDAC with matched liver metastases
to better understand the roles of ITH and TME in PDAC evolution at the
single-cell resolution. Our study may provide an exciting entry point
for developing novel therapeutic strategies taking into account the ITH
and TME of PDAC.
Supplementary Information
[174]Supplementary Material 1.^ (6.9MB, pdf)
Acknowledgements