In total, 1,109 PanNET patients from the multicentric database were included in this study (Supplement Table 1). The proportions of G1, G2 and G3 were 43.1%, 51.0%, and 5.9%, respectively. A total of 30.0% of patients had distant metastasis before surgery. The median survival time was 144.0 months (5-year survival rate, 82.7% [74.8%~88.3%]; 10-year survival rate, 68.9% [52.1%~79.6%]).
To explore risk factors for PanNET, multivariate Cox analysis was performed. After excluding collinear factors, the AJCC grading was dependent on metastasis, and multivariate Cox proportional hazards regression analysis identified G grade and AJCC M stage as independent risk factors for OS (Table 1). However, cross-tabulation analysis showed only a weak correlation between metastasis and G grade (Spearman rho = 0.09, P = 0.012, Supplement Table 2). Furthermore, there was a high proportion of G1 (31.88%) and G2 (59.42%) patients with metastasis.
Table 1
The univariate and multivariate analysis of overall survival in 1109 PanNET patients. The univariate analysis results showed that AJCC grading, Pathology grading, AJCC T stage, AJCC N stage, AJCC M stage, lymph node metastasis, tumor location could be considered as influence factors in panNEN survival. In the further multivariate analysis for overall survival, AJCC T stage, AJCC M stage, Pathology grading remained as independent risk indicators for tumor OS.
Variable
|
Univariate analysis
|
Multivariate analysis
|
HR(95%CI)
|
P
|
HR(95%CI)
|
P
|
Age ≥ 60
|
1.206(0.902 ~ 1.611)
|
0.2063
|
|
|
Female
|
0.757(0.572 ~ 1.001)
|
0.0512
|
0.850 (0.638 ~ 1.133)
|
0.2678
|
G grade
|
|
|
|
|
G1
|
Reference
|
|
Reference
|
|
G2
|
1.611(1.171 ~ 2.217)
|
0.0034
|
1.134 (0.816 ~ 1.575)
|
0.4531
|
G3
|
7.600(4.859 ~ 11.889)
|
< .0001
|
5.293 (3.276 ~ 8.553)
|
< .0001
|
T Stage
|
|
|
|
|
T1
|
Reference
|
|
Reference
|
|
T2
|
3.708(2.088 ~ 6.584)
|
< .0001
|
3.138 (1.756 ~ 5.606)
|
0.0001
|
T3
|
7.359(4.192 ~ 12.917)
|
< .0001
|
3.519 (1.937 ~ 6.394)
|
< .0001
|
T4
|
6.566(3.386 ~ 12.729)
|
< .0001
|
3.519 (1.937 ~ 6.394)
|
0.0001
|
M1 Stage
|
8.085(5.801 ~ 11.269)
|
< .0001
|
6.071 (4.251 ~ 8.670)
|
< .0001
|
N stage
|
|
|
|
|
N0
|
Reference
|
|
Reference
|
|
N1
|
3.740(2.825 ~ 4.950)
|
< .0001
|
2.079 (1.523 ~ 2.839)
|
< .0001
|
AJCC stage (8th edition)
|
|
|
|
|
I
|
Reference
|
|
|
|
II
|
3.328(1.806 ~ 6.133)
|
0.0001
|
|
|
III
|
7.304(3.940 ~ 13.540)
|
< .0001
|
|
|
V
|
30.637(16.158 ~ 58.089)
|
< .0001
|
|
|
Primary Site
|
|
|
|
|
body/tail
|
Reference
|
|
Reference
|
|
head
|
1.492(1.128 ~ 1.972)
|
0.005
|
1.482 (1.101 ~ 1.996)
|
0.0095
|
other
|
0.688(0.216 ~ 2.192)
|
0.5276
|
0.642 (0.197 ~ 2.088)
|
0.4614
|
As both G grade and metastasis are important prognostic factors for PanNETs, we combined them and explored the ability to predict survival. Thereafter, six subgroups according to the three G stage and two M stage (G 1–3/M 0–1) were entered into the multivariate model along with the main effects (Supplementary Table 3). Based on the median survival time and HR, a new G/M staging classification including three groups was derived (Fig. 1A-B): G/M-I, including G1M0 and G2M0, with a median survival time of 179.2 months (HR < 2); G/M-II, including G1M1, G2M1 and G3M0 (median OS survival: 81.0 months; HR: 7–12); and G/M-III, including G3M1 (median OS: 12.0 months; HR = 59.3). The survival time of G1- and G2-grade patients with metastasis was similar to that of G3-grade patients without metastasis (Supplementary Table 4).
The discrimination ability of the new grouping was evaluated with the Harrell C-index. The apparent C-index was 0.76 (95% CI: 0.73–0.80), indicating better discrimination ability than G grade (0.61 [95% CI: 0.57–0.65]) or metastasis (0.66 [95% CI: 0.62–0.70]) alone. The predictive accuracy using time-dependent ROC analysis was relatively high at 0.747 (95% CI: 0.677–0.816) at 3 years and 0.726 (95% CI: 0.661–0.791) at 5 years (Fig. 1C). These findings suggest that G grade and metastasis together can more accurately predict OS. G grade may also be evaluated by Ki-67. However, there are no suitable markers to evaluate the metastatic ability of PanNETs, and screening suitable markers is thus necessary.
Single-cell RNA-seq highlights PanNET characteristics of metastasis do not correspond with G grade.
To explore the existence of high metastasis potential subgroups in PanNET, we generated single-cell RNA-seq profiles from 4 PanNET tumor samples and 3 normal pancreas (NP) samples without any treatment performed. After initial quality control, we acquired single-cell transcriptomes for a total of 39,374 cells from PanNET samples and 17,711 cells from control pancreas samples. To define the cellular composition in the tumor microenvironment and subpopulation structure of PanNET, we performed graph-based clustering analysis on these cells and identified 14 distinct cell types based on known marker genes or specific genes expressed by each cell type, including EC, cancer cells, ductal cells, stellate cells, acinar cells, fibroblasts, DCs, macrophages, monocytes, T cells, endocrine cells, Langerhans cells, B cells, and mast cells (Fig. 2D, Supplementary Fig. S1A).
We observed great differences in the cell distribution and cell fraction for each cell type for these samples from 4 patients (Fig. 2E, 2F). In patients 1 and 4, the proportion of tumor cells to all cells was 94.7% and 96.3%, respectively; the proportion was 8.8% and 9.8%, respectively, in patients 2 and 3,with a high concentration of endothelial and stellate cells (Fig. 2F).
Large-scale chromosomal copy number variation (CNV) analysis for 4 patients showed that PanNET1 and PanNET4 exhibited markedly higher CNV levels than PanNET2 and PanNET3 (Fig. 2G, supplement Fig S1B). Next, the intertumoral cellular heterogeneity score (ITH score) showed that PanNET2 and PanNET3 had higher intertumoral cellular heterogeneity (Fig. 2H). With regard to entropy, PanNET4 had the lowest score, indicating a poor differentiation level (Fig. 2I). Although PanNET1, PanNET2 and PanNET3 had similar scores, G grades differed (Fig. 2I). Transcription factors also varied in the four tumors (Supplementary Fig S1C and S1D). These results indicate significant differences between the 4 tumor samples.
GSVA revealed activated pathways related to invasion and metastasis in PanNET2 and PanNET3. PanNET4 exhibited a weak invasive capability (Fig. 2J). Interestingly, PanNET1-3 were relatively aggressive, whereas PanNET4 has a weak invasive capability, similar to an indolent tumor. GO analysis confirmed activation of the focal adhesion pathway and cell-substrate junction pathway, which are both closely related to the invasion process, in PanNET2 and PanNET3 (Supplementary Fig S1E-H).
These observations indicate that PanNET characteristics of metastasis do not correspond with G grade.
Single-cell RNA-seq identifies high metastatic capacity subgroups in PanNETs.
Cancer cells were further divided into 16 subgroups based on t-SNE analysis (Fig. 2A, 2B). By comparing gene expression levels, we found each subgroup to express a specific gene set that can be used for distinguishing these subgroups (Supplementary Fig S2A). Notably, 4 PanNETs had a different subgroup composition (Fig. 2C).
Large-scale chromosomal CNVs were found in all subgroups, and different subgroups exhibited different CNV levels (Fig. 2D, Supplementary Fig. S2B). Subgroups 11 and 13 exhibited remarkably higher ITH scores than the other subgroups (Fig. 2E), and subtypes 9, 11 and 13 exhibited remarkably higher levels of entropy (Fig. 2F). Transcription factors also varied among the 16 subgroups (Supplementary Fig S2C).
Additionally, GSVA showed activation of pathways related to invasion and metastasis in subgroups 9, 11 and 13 (Fig. 2G). Subgroup 11 was found to be the main component of PanNET2, and subgroups 9 and 13 were the main components of PanNET3, which is consistent with the results of GSVA analysis for these two tumors. We then performed functional enrichment analysis for each subgroup and found subgroups 9, 11 and 13 to be enriched in migration-related terms, such as response to focal adhesion, emphasizing their potential functions in migration and metastasis (Supplementary Fig S2F and S2G).
In conclusion, the subgroup composition of tumors varied, even for the same G-grade tumors. We identified 3 subgroups with high metastatic potential, and tumors with a high proportion of these high metastatic potential subsets have a stronger ability to metastasize.
PIN1 is highly expressed in PanNETs and associated with metastasis.
To screen metastasis markers in PanNETs, we selected 32 surgically resected primary PanNET tumor tissues and 10 nontumor pancreas tissues for proteomics analysis (Supplemental Fig S3A, Supplemental Table 5).
A total of 9799 protein groups were quantified from 32 PanNET tumor samples and 10 adjacent nontumor (NT) samples by the DIA proteomic approach (Fig. 3A). In total, 1628 proteins were significantly altered between PanNET and NT, of which 882 and 746 were upregulated and downregulated, respectively, in PanNET tumor tissues (Supplemental Fig S3B). Heatmap and PCA of the significantly regulated proteins showed that the majority of samples displayed similar protein expression patterns and could be grouped into the same biological conditions (Fig. 3B, supplement Fig S3A). GO analysis showed that proteins involved in metastatic processes, such as cell adhesion, extracellular matrix organization and collagen fibril organization, were enriched for these upregulated and downregulated proteins. According to KEGG analysis, the upregulated and downregulated proteins are involved in pathways such as focal adhesion and ECM-receptor interaction (Fig. 3C, supplement Fig S3C).
To screen for proteins regulating metastasis in PanNETs, 10,154 protein groups were quantified from the nonmetastasis group and distant metastasis group (Supplemental Fig S3D). A total of 225 proteins were significantly altered, of which 108 were upregulated and 117 downregulated in metastatic tissues compared to nonmetastatic tissues (Fig. 3D, supplement Fig S3D-E). GO analysis revealed enrichment of biological processes, such as invasion and cell adhesion, among the proteins upregulated and downregulated in the metastasis group. KEGG analysis showed pathways such as morphine addiction and focal adhesion to be enriched among the upregulated proteins (Fig. 3E, supplement Fig S3F).
Next, we found 18 proteins that were both highly expressed in PanNET tumors and in tumors with metastasis (Fig. 3F, supplement Fig S3G-H). Expression of these 18 genes in single-cell RNA-seq data showed PIN1, POSTN, SEMA4F, ASPN and KCDT12 to be relatively highly expressed in tumor cells (Supplemental Fig S3I). PIN1, POSTN, SEMA4F, ASPN and KCDT12 are related to tumor invasion and metastasis according to a published article and were selected for further study.
Then, immunohistochemical assays of the five selected proteins (Supplemental Fig S4A) were carried out using 152 PanNET specimens (Supplemental Table 6). Among the five proteins, only PIN1 and POSTN expression had significant prognostic value (Fig. 4A-D, supplement Fig S4B-G). Furthermore, PIN1 and POSTN showed stronger expression in 29 tumors from patients who died during follow-up than in patients who survived (Fig. 4E). PIN1 and POSTN showed stronger expression in 44 tumors with recurrence than in 104 tumors without recurrence (Fig. 4F).
To further examine the importance of the independent and select metastasis-specific variables, we applied a nonparametric approach for competing risks using random survival forests and that PIN1 had the highest variable importance (100.0) associated with metastasis (Fig. 4G): the higher the expression of PIN1 was, the more distant the metastasis was (Spearman correlation coefficient = 0.30, P < 0.0001). For survival-specific variables, we applied a nonparametric approach for competing risks using random survival forests and found that PIN1 and G grade significantly affected OS (VIMP > 0.05, Fig. 4H). In addition, POSTN and PIN1 correlated linearly and directly (P interaction = 0.025). To exclude multicollinearity, only PIN1, which is more closely related to metastasis, was assessed in multivariate analysis, which showed that PIN1 and G stage affected OS (Table 2). Next, we analyzed the correlation between Ki-67 and PIN1 expression in PanNET, and the results showed no linear correlation (Spearman correlation coefficient = 0.14, P = 0.778). (Supplemental Fig S4H-J)
These data suggest that PIN1 is highly expressed in PanNETs and may serve as a marker of the metastatic capability of PanNETs.
Lentivirus and sulfopin1 were used in QGP-1 cell lines to affect the activity of PIN1, evaluating whether the proliferation and migration of PanNET were affected by PIN1. RNA-seq and GSEA analysis were conducted and showed that the use of sulfopin1 was related to the down-regulation of EMT pathway in QGP1, which was the most inhibited pathway (Fig. 5A, B, D). and overexpression of PIN1 was related to the up-regulation of EMT pathway in QGP1 (Fig. 5C). Further, wound healing assay suggested that using PIN1 inhibitors inhibited the migration of QGP1 (Fig. 5E). Overexpression of PIN1 up-regulated Vimentin, down-regulated E-cadherin and promoted S675 phosphorylation of CTNNB1, suggesting the activation of EMT pathway (Fig. 5F). On the contrary, knocking down PIN1 and using sulfopin1 down-regulated Vimentin, up-regulated E-cadherin and inhibited S675 phosphorylation of CTNNB1, suggesting inhibition of EMT pathway (Fig. 5G, H).
In order to clarify the molecular mechanism of PIN1 regulating the proliferation and migration of PanNET cells, IP-MS analysis of PIN1 was conducted and showed that 97 proteins interacted with PIN1 (Fig. 5I). Using the data of proteome and phosphorylation omics from the CPTAC database for bioinformatics analysis, the correlation between proteins or phosphorylation sites and PIN1 expression were analyzed, and six proteins (SNRNP200, PUM1, HSP90AB1, ANK3, RPLP2, LMNA) were found to have potential interactions with PIN1 (Fig. 5J). The interaction between PIN1 and LMNA was verified by exogenous co-immunoprecipitation (Fig. 5K, L). Since the activity of proline isomerase of PIN1 depends on the phosphorylation of serine or threonine on the T/S-P motif near the N-terminal of the target protein, the conservation of LMNA amino acid sequences among species was further verified (Fig. 5M). Serine and threonine at T/S-P motif on LMNA were mutated to alanine one by one, and the co-immunoprecipitation analysis showed that the binding between PIN1 and LMNA was significantly inhibited by LMNA-S22A mutation (Fig. 5N).
To explore the effect of PIN1 on clinical treatment decisions, we established a survival signature that can classify 152 PanNET patients into homogeneous subpopulations according to OS. After pruning the decision trees using the postpruning method, 4 terminal nodes (subpopulations) representing a survival signature were identified (Fig. 6A). Patients harboring the PIN1lowG1 signature (Node 4) had a 5-year OS rate of 96.8%; in contrast, patients harboring the G3 signature (Node 7) had the worst 5-year OS rate (0%). A survival index based on a multivariate Cox proportional hazards model of independent prognostic variables and OS was calculated for each individual in the 4 subpopulations (Fig. 6B).
The waterfall plot of survival indices within the subpopulations indicated relevance between the survival signature and the survival index (Fig. 6C). The discrimination ability of the former was evaluated with the Harrell C-index. The apparent C-index was 0.921 (95% CI: 0.859–0.983) and the bias-corrected C-index using a bootstrapping procedure (1000 resamples) was 0.898 (95% CI: 0.773–0.984), indicating a better discrimination ability than G grade (0.77 [95% CI: 0.68–0.86]) or metastasis (0.61[95% CI: 0.52–0.69]) alone.The predictive accuracy of the survival signature using time-dependent ROC analysis was relatively high at was 0.933 (95% CI: 0.876–0.990) at 3 years and 0.909 (95% CI: 0.837–0.981) at 5 years (Fig. 6D). Calibration curves also showed good agreement between the predicted and observed OS (Fig. 6E).