
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)02061-3
10.1016/j.isci.2024.110836
110836
Article
Single-cell sequencing reveals the heterogeneity of pancreatic neuroendocrine tumors under genomic instability and histological grading
Ye Zeng 1234515
Zhou Yan 615
Hu Yuheng 715
Li Qiang 815
Xu Zijin 915
Lou Xin 12345
Zhang Wuhu 12345
Zhu Di 10
Xie Cao 11
Zhou Qingtong 1012
Gao Jing 13
Zhou Hu 13
Yang Dehua 612
Qin Yi 12345
Xu Xiaowu 12345
Chen Jie chenjie@fudanpci.org
2∗
Ji Shunrong jishunrong@fudanpci.org
1234516∗∗
Wang Ming-Wei mwwang@simm.ac.cn
101214∗∗∗
Yu Xianjun yuxianjun@fudanpci.org
12345∗∗∗∗
1 Department of Pancreatic Surgery, Fudan University Shanghai Cancer Center, Shanghai 200032, China
2 Center for Neuroendocrine Tumors, Fudan University Shanghai Cancer Center, Shanghai 200032, China
3 Department of Oncology, Shanghai Medical College, Fudan University, Shanghai 200032, China
4 Shanghai Pancreatic Cancer Institute, Shanghai 200032, China
5 Pancreatic Cancer Institute, Fudan University, Shanghai 200032, China
6 The National Center for Drug Screening, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China
7 Department of Hepatobiliary and Pancreatic Surgery, Tenth People’s Hospital of Tongji University, School of Medicine, Tongji University, Shanghai, China
8 Department of General, Visceral, and Transplant Surgery, Ludwig-Maximilians-University Munich, Marchioninistr. 15 81377 Munich, Germany
9 General Surgery Department, Qingpu Branch of Zhongshan Hospital, Shanghai 200032, China
10 Department of Pharmacology, School of Basic Medical Sciences, Fudan University, Shanghai 200032, China
11 School of Pharmacy, Fudan University, Shanghai 201203, China
12 Research Center for Deepsea Bioresources, Sanya, Hainan 572025, China
13 The CAS Key Laboratory of Receptor Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China
14 Department of Chemistry, School of Science, The University of Tokyo, Tokyo 113-0033, Japan
∗ Corresponding author chenjie@fudanpci.org
∗∗ Corresponding author jishunrong@fudanpci.org
∗∗∗ Corresponding author mwwang@simm.ac.cn
∗∗∗∗ Corresponding author yuxianjun@fudanpci.org
15 These authors contributed equally

16 Lead contact

27 8 2024
20 9 2024
27 8 2024
27 9 1108361 9 2023
10 1 2024
23 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Histological grading is the key factors affecting the prognosis and instructive in guiding treatment and assessing recurrence in non-functional pancreatic neuroendocrine tumor (NF-Pan-NET). Approximately one-third of patients without copy number variation (CNV) alteration and the prognosis of these patients are better than that of patients with CNV alteration. However, the difference between CNV and histological grading is unclear. Here, we analyzed the heterogeneity of tumor cells according to two classification criteria, genomic instability (including CNV alteration and tumor mutation burden) and histological grading. We revealed that the activated core pathways of tumor cells were significantly different under different histological grading’s and genomic instability patterns. We also found that tip cells, lymphatic endothelial cells, macrophages, CD1A + dendritic cell, Treg, MAIT, ILC, and CAFs might participate in the process of hepatic metastases, which will facilitate the understanding of the patterns to decode the malignant potential and of NF-Pan-NET.

Graphical abstract

Highlights

• Pan-NET tumor cells show marked heterogeneity in copy number variation

• mTOR signaling pathway is activated in CNV-altered tumor cells

• CNV neutral tumor cells expressed the characteristics of normal neuroendocrine cells

• Metastatic tumor cells have a unique gene expression profile

Microenvironment; Cancer; Omics; Transcriptomics

Subject areas

Microenvironment
Cancer
Omics
Transcriptomics
Published: August 27, 2024
==== Body
pmcIntroduction

Pancreatic neuroendocrine neoplasm (Pan-NEN) ranks the second most common malignancies of the pancreas, which originate from the neuroendocrine precursor cells or neuroendocrine transdifferentiation of epithelial cell.1,2 According to the data from the Surveillance, Epidemiology, and End Results (SEER) program, in the context of a plateau in the overall incidence of all malignant neoplasms, the incidence of NEN increased by a staggering 7 times between 1973 and 2012 in the United States.3 Pan-NEN account for 30% of gastroenteropancreatic neuroendocrine neoplasm (GEP-NEN) and 1–2% of all pancreatic malignancies.4 Compared with other GEP-NEN, Pan-NEN has a higher rate of liver metastasis, up to 20%–64%, and diffuse liver metastasis was the most common.5,6 It is characterized by high heterogeneity, and can manifest slow inert growth, invasive growth or even early metastasis. These characteristics may vary with the progression of the disease. According to cell differentiation, well-differentiated was defined as Pan-NET and poorly differentiated was defined as pancreatic neuroendocrine carcinoma (Pan-NEC). Pan-NET was divided into functional and non-functional Pan-NET (NF-Pan-NET) according to whether it secreted hormones and caused clinical symptoms. Further, based on mitotic image and ki-67 proliferation index, Pan-NET were divided into three grades of G1, G2, and G3. Histological grading is instructive in guiding treatment, assessing recurrence, and estimating prognosis. Generally, G1 indicates the best prognosis and G3 implies the worst prognosis.4 However, once G1 and G2 Pan-NET are complicated by liver metastases, the prognosis is extremely poor, the median survival is 24 months7 Therefore, decoding the liver recurrence and metastasis of Pan-NEN can radically improve the prognosis. In addition, the histological grading is not static. In some cases, the grade increases with the progression of the disease, especially in the metastases.8 Increased grade tends to predict a worse prognosis.9 Besides, 90% of Pan-NET are grade G1 or G2, leading to an unpredictable clinical course that varies from indolent to highly malignant. Therefore, there is an urgent need to find better patterns to decode the malignant potential of Pan-NET.

One previous study recruited whole-genome sequencing (WGS) of 98 Pan-NET, in which the authors revealed four discrete groups of patients based on arm length copy number patterns, including: (1) recurrent pattern of whole chromosomal loss (RPCL); (2) limited copy number events, many of which were losses affecting chromosome 11; (3) polyploidy; and (4) aneuploidy.10 Increased genomic instability such as increased tumor mutation burden (TMB) and copy number variations (CNV) is the basis of tumor development. Most of tumors show high levels of mutations or high levels of CNVs, and some belongs to an intermediate level of both.11 However, according to the data of these 98 Pan-NET, more than a third of the population, possessed limited copy number events and lowest mutation burden.10 This suggests that these patients have low genomic instability. Another study recruited whole-genome/whole-exome sequencing (WGS/whole exome sequencing [WES]) of 211 Pan-NET also revealed novel categories of Pan-NET based on CNV patterns: amplification, copy neutral, and deletion. It was found that non-functional pancreatic neuroendocrine tumors (NF-Pan-NET) with CNV alterations (amplification and deletion) had an elevated risk of relapse and worse survival. Same as the previous study, the group of copy neutral had the lowest TMB. In addition, the distribution of histological grading also varied from the CNV patterns: in the CNV neutral group, G1 accounted for more than 50%, while in the CNV amplification/deletion groups, G2 dominated.12 This suggests that the combination of histological grading and CNV alterations classification system may have important significance for the treatment and prognosis assessment of NF-Pan-NET.

Single-cell sequencing has unique superiorities in interpreting tumor heterogeneity. Chromosome instability, found in almost every tumor, is the internal impetus for tumor evolution, endowing tumors with multifarious clones. This internal impetus coupled with tumor microenvironment provokes high heterogeneity intra or inter-tumors.13,14 However, the commonly used bulk sequencing, which profiles a holistic genetic landscape of many cells, shows obviously deficiencies in disclosing tumor heterogeneity. In this study, we recruited 17 NF-Pan-NET patients and analyzed their single-cell sequencing profiles to map associated genetic landscape and cell type. We further explored the heterogeneity in NF-Pan-NET with different histological grading and genomic instability patterns, in addition to identifying the characteristics of patients with hepatic metastases (PHM).

Results

Single-cell profiling of the tumor ecosystem in NF-Pan-NET

To explore the tumor ecosystem in NF-Pan-NET, we performed scRNA-seq using the Chromium platform (10× Genomics) from 15 NF-Pan-NET tumor specimens (primary tumor tissues, PT), 3 control pancreas (adjacent normal tissues, NT) and 2 hepatic metastases (HM) without any treatment. Besides, further WES, bulk RNA-sequencing and proteomics were performed in six or seven of these patients (Figure 1A). The demographic and pathological data of the 17 randomly selected NF-Pan-NET cases were shown in Figure 1A. To better identify tumor cells, we also included scRNA-seq data of 8 control pancreases from the study by Peng et al. (Numbers N1, N4-9, and N11 in the original study).15 These data were used to characterize tumor heterogeneity and hepatic metastases in NF-Pan-NET (Figure 1B). After the initial quality control, we acquired scRNA-seq data of 123 106 cells in total. Uniform manifold approximation and projection (UMAP) algorithm was performed to reduce nonlinear dimensionality and visualize these data. These cells were clustered into 47 clusters corresponding to 12 cell types, including B cells, T cells, smooth muscle cells (SMC), plasma cells, mast cells, macrophages, fibroblasts, endothelial, endocrine cells, ductal cells, cycling cells, and acinar cells (Figure 1C). Marker genes of each cluster were calculated by FindAllMarkers function with wilcox rank-sum test algorithm (Figure 1D). Furthermore, CIBERSORTx was used to analyze the proportion of these 12 cell types in 83 bulk RNA-seq samples from Alvarez’ study.16 Comparing with the proportion in our scRNA-seq data, we found that for half of the cell types in the two datasets, there was no significant difference in their proportion. However, in our data, the proportion of endocrine cells was significantly reduced (although still the most abundant), while the proportion of some immune cells was significantly increased (B cells, macrophages, mast cells, and T cells) (Figure 1E). In addition, the proportion of endocrine cells, B cells, macrophages, T cells, and mast cells were increased in primary tumor tissues while ductal cells, acinar cells, endothelia cells, and fibroblast showed an obviously opposite trend in our data (Figure 1E). Intriguingly, there was no significant difference in the distribution of these cells between the G1 and G2 (Figure 1F). In general, tumor tissue and normal tissue have significantly different dominant cell populations (Figure 1G). These data initially delineated the single-cell profiling of the tumor ecosystem in NF-Pan-NET. Notably, there was no endocrine cell found in T4 and T6 samples. Therefore, these two specimens will be excluded when analyzing the heterogeneity of tumor cells.Figure 1 Single-cell profiling of the tumor ecosystem in NF-Pan-NET

(A) The demographic and pathological data of the 17 randomly selected NF-Pan-NET cases in our center, single-cell RNA-sequence, and multiple-omics (WES, RNA-sequencing and proteomics) of bulk tissues were conducted in some of cases.

(B) Schematic diagram showed the design and purpose of this study.

(C) The cells were visualized by UAMP.

(D) Marker genes of each cluster were shown in the heatmap.

(E) The proportion of selected cell types in bulk RNA-seq was calculated by CIBERSORTx (left) and the proportion of selected cell types in our enrolled scRNA-seq data. The data of scRNA-seq are represented as mean ± SEM. The statistical analysis in the left showed the comparation of cell proportion between bulk RNA-seq and scRNA-seq data. The statistical analysis in the right showed the comparation of cell proportion in scRNA-seq data between group NT and PT.

(F) Proportion of cell types in different groups.

(G) Heatmap showed the proportion of cell types in each specimen.

Tumor cells show marked heterogeneity in copy number variation

To define malignant cells, we calculated large-scale chromosomal copy number variation (CNV) in endocrine cells by InferCNV. We found that endocrine cells in tumor tissues exhibit significant differences in CNV patterns. In T1 (3.94%), T11 (4.98%), and T13 (23.33%) the malignant cells identified by inferCNV were only a small part of the endocrine cells in the tissues. In H1 (36.85%), H2 (43.18%), and T15 (49.89%), which were from the same patient, and in T3 (40.30%) and T12 (41.68%) the malignant cells identified by inferCNV accounted for a medium portion of the endocrine cells in the tissues. While in T2 (96.27%), T5 (91.76%), T7 (99.54%), T8 (83.12%), T9 (92.18%), T10 (99.33%), and T14 (56.35%) the malignant cells identified by inferCNV accounted for the vast majority of the endocrine cells in the tissue (Figures 2A, 2B, and S1A). To further confirm the CNV differences, we performed WES in T2, T3, T6, T10, T11, and T14 bulk tissues, and we found that T3 and T11 had neutral CNV, T6 had a deleted CNV, whereas T2, T10, and T14 had amplified CNV (Figure 2C). Intriguingly, T3 and T11 presented neutral CNV in bulk tissues and had lower proportions of inferCNV-high neuroendocrine cells in scRNA-sequence. Correspondingly, T2, T10, and T14 presented amplified CNV in bulk tissues and had higher proportions of inferCNV-high neuroendocrine cells in scRNA-sequence. These indicated that there was a strong correlation between CNV in the WES data of bulk tissues and the proportion of inferCNV-high neuroendocrine cells in scRNA-sequence. We further calculated the TMB using the WES data. We found that T6, T2, T10, and T14 had higher TMB than that of T3 and T11 (Figure 2D). Both CNV and TMB are important elements of genomic instability, which play an important role in the development of cancers. Generally, there is a striking inverse relationship between CNV and TMB at the extremes of genomic instability, particularly in highly altered tumors.11 However, in NF-Pan-NET, the patients with CNV alteration tended to have higher TMB and worse prognosis,10,12 suggesting that there is great heterogeneity in the pathogenesis of Pan-NET. However, single-cell profiling in the context of genomic instability pattern of NF-Pan-NET is unknown. The proportion of normal pancreatic endocrine cells in the pancreatic tissue is very low (Figure 1F). Therefore, the inferCNV-low endocrine cells still should be tumor cells. To confirm it, we analyzed the signal pathway of endocrine cells in different specimens. We found that the endocrine cells in tumor tissues had activated signal pathways, including DNA repair, ROS, mTOR signaling, and so on. However, endocrine cells in normal pancreas had activated signal pathways associated with endocrine and even exocrine functions of the pancreas, such as insulin secretion, bile secretion, fat digestion and absorption, protein digestion and absorption, pancreatic secretion, and hallmark of pancreas beta cells (Figures 2E and S1B). To further verify this, we analyzed the data of transcriptomics and proteomics in bulk-tissues. We calculated the tumor purity through the transcriptome data and found that the selected sample had a high tumor purity (Figure S1C). Our results showed that tumor tissues exhibited reduced normal pancreatic function and enhanced tumor-related pathways regardless of genomic stability (Figures S1D and S1E). These results suggested that the identified endocrine cells in tumor tissues by scRNA-seq were tumor cells. In addition, tumor cell showed marked heterogeneity in genome instability and inferCNV had limitation in identifying tumor cells in scRNA-seq data of NF-Pan-NET. We defined the patients with high CNV alteration (CNV deletion and CNV amplification) and high TMB as high genome instability group (GIH). Correspondingly, patients with low CNV alteration (CNV neutral) and low TMB as low genome instability group (GIL). We defined T1, T3, T11, and T13, which possessed a low proportion of inferCNV-high neuroendocrine cells, as the group of low genome instability (GIL). Correspondingly, T2, T5, T7, T8, T9, T10, T12, and T14, which possessed a high proportion of inferCNV-high neuroendocrine cells, as the group of high genome instability (GIH). Notably, despite the malignant cells identified by inferCNV accounted for a medium portion of the endocrine cells in in T3 (40.30%) and T12 (41.10%), T3 was classified in the GIL group due to the bulk-sequence indicating that it presented neutral CNV. T12 was classified in the GIH group due to its high CNV score, which was similar to T2, T7, T8, T9, and T10. The proportion of malignant cells identified by inferCNV was 56.35% in T14. T14 was classified in the GIH group due to its bulk-sequence indicates that T14 presented amplified CNV. Since H1 (36.85%), H2 (43.18%), and T15 (49.89%) were from the same patient and in these specimens, we defined these special tissues as group of NF-Pan-NET with hepatic metastases (PHM) (Figure 2A).Figure 2 Tumor cells show marked heterogeneity in copy number variation

(A) CNV score of all endocrine cells in each sample. The ordinate represents the frequency.

(B) CNV score of selected CNV-high endocrine cells in each sample. The ordinate represents the frequency. Percentage represents the proportion of inferCNV-high neuroendocrine cells to total neuroendocrine cells.

(C) Copy number variations in the WES data.

(D) TMB were analyzed in the WES data.

(E) Signal pathway of endocrine cells in different specimens.

Revealing the heterogeneity of tumor cells under the pattern of genomic instability and histological grading

To further unscramble the heterogeneity of tumor cells, we tried to explore the rules of tumor cell subtype analysis. We depicted all endocrine cells in the tumor tissue and normal pancreas by t-distributed stochastic neighbor embedding (tSNE) map and found 18 sub-clusters of neuroendocrine cells (Figure 3A). Histological grading, based on the mitotic image and ki-67 proliferation index, is instructive in guiding treatment, assessing recurrence, and estimating prognosis in NF-Pan-NET.4 However, in the tSNE map of endocrine cells, we found that there was no obvious distribution pattern of tumor cells between G1 and G2 (Figure 3B). We further found that the distribution of tumor cells in this particular patient specimen with hepatic metastases was significantly different from that in other patients. Besides, in this patient, there was also no clear demarcation between the hepatic metastases and the primary tumor (Figure 3C). Thus, it was reasonable to assign this patient’s tumor cells to a separate group (PHM). However, when inferCNV values were used as a classification basis, tumor cells were clearly differentiated (Figure 3D). To further demonstrate the role of CNV, we calculated the proportion of tumor cells with high CNV in each cluster. We found that cluster 4, 5, 8, 9, and 12 were mainly composed of CNVhigh tumor cells. Cluster 7, 10, 11, 13, 14, 15, and 16 were mainly composed of CNVlow tumor cells. However, in cluster 0, 1, 2, 3, 6, and 17 CNVhigh tumor cells accounted for nearly half (Figure 3E). We further analyzed the subcluster of tumor cells according to the source of tumor cells and the group of patients. We found that cluster 4, 5, 8, 9, and 12 were mainly composed of tumor cells from patients with G1 grade. Cluster 7, 10, 11, 13, 14, 15, and 16 were mainly composed of tumor cells from patients with G2 grade (Figure 3F). Intriguingly, if the sample was classified according to genomic instability, cluster 4, 5, 8, 9, and 12 were almost composed of tumor cells from patients in GIH group and cluster 7, 10, 11, 13, 14, 15, and 16 were almost composed of tumor cells from patients in GIL group. Cluster 0, 1, 2, 3, 6, and 17 were almost composed of tumor cells from the PHM (Figure 3G). This indicated that genomic instability was more suitable for classifying tumor cells than histological grading. In addition, the proportion of tumors from group GIH in the subcluster was positively correlated with the proportion of tumors from patients with G1 grade (Figure S2A). Similarly, the proportion of tumors from group GIL in the subcluster was positively correlated with the proportion of tumors from patients with G2 grade (Figure S2B). For patients with liver metastases, there was no significant relationship between the proportion of CNVhigh tumor cells in the subcluster between the proportion of tumors from primary or metastatic site (Figure S2C). However, the proportion of CNVhigh tumor cells in the subcluster was positively correlated with the proportion of tumors from GIH group and patients with G1 grade and was negatively correlated with the proportion of tumors from GIL group and patients with G2 grade (Figure S2D). Further analysis by CytoTRACE showed that cluster 0, 1, 2, 3, 6, and 17, which were mainly from group PHM, exhibited the similar differentiation state, followed by cluster 4, 5, 8, 9, and 12, which were mainly from group GIH. Cluster 7, 10, 11, 13, 14, 15, and 16, which were mainly from group GIL, exhibited the state of terminal differentiation, which was similar to normal endocrine cells (Figure 3H). Gene-module analysis showed that cluster 4, 5, 8, 9, and 12, which were mainly from group GIH, had similar gene expression mode, such as mTOR signaling pathway; cluster 0, 1, 2, 3, 6, and 17, which were mainly from group PHM, had similar gene expression mode, such as TGF-β signaling pathway; sub-clusters 11, 14, 15, and 16, which mainly from group GIL, had similar gene expression mode such as Notch signaling pathway; cluster 7, 10, and 13, and, which were mainly from group GIL and normal neuroendocrine cells, had similar gene expression mode such as Wnt signaling pathway and insulin secretion (Figure 3I). We further analyzed the degree of pathway activation using quantitative set analysis for gene expression (QuSAGE) and found that tumor cell from group PHM exhibited activated hedgehog signaling, TGF-β signaling, glycolysis, hypoxia, and Wnt signaling pathway; tumor cell from group GIH exhibited activated DNA repair, ROS, E2F, and growth-factor signaling; while tumor cell from group GIL exhibited activated pancreatic beta cells and pentose phosphate pathway (Figure 3J). Tumor cells from patients with G1 grade exhibited activated Wnt-β-catenin signaling pathway; while tumor cells from patients with G2 grade exhibited activated ROS signaling (Figure 3K).Figure 3 Revealing the heterogeneity of tumor cells under the pattern of genomic instability and histological grading

(A) The subcluster of endocrine cells was presented on the tSNE map.

(B–D) The source of endocrine cells was presented on the tSNE map. PT-G1, endocrine cells from the patients with histological grading G1. PT-G2, endocrine cells from patients with histological grading G2. HM, hepatic metastases. NT-a, endocrine cells of normal pancreas from patients in our center. NT-b, endocrine cells of normal pancreas from Peng’s study. PT non-metastatic, endocrine cells from patients without hepatic metastases. PT metastatic, endocrine cells of primary lesion from patients with hepatic metastases.

(E) The proportion of CNVhigh tumor cells, CNVlow tumor cells, and normal endocrine cells in each cluster.

(F) The proportion of tumor cells from patients with G1 grade, G2 grade, hepatic metastases, and normal adjacent tissues in each cluster.

(G) The proportion of tumor cells from group GIL, GIH, PHM, and normal adjacent tissues in each cluster.

(H) Differentiated degree of endocrine cells in each cluster was analyzed by CytoTRACE.

(I) Gene-module analysis of endocrine cells in each cluster.

(J and K) The activation of KEGG pathway in endocrine cells of different group was analyzed by QuSAGE.

In order to further explore the relationship between genomic instability and histological grading, we perform WES on 77 NF-Pan-NET tumor specimens and analyzed the CNV alteration (Figure S2E). We divided the patients into 3 group based on the overall CNV and the presence of hepatic metastases (PHM, GIH, and GIL). Indeed, in group PHM, 5 out of 11 patients showed low genome instability and 6 out of 11 patients showed high genome instability (Figure S2E). In this cohort, patients in group GIL or patients with G1 grade showed significantly longer relapse-free survival (RFS) compared with patients in group GIH or patients with G2 grade (Figures S2F and S2G). Since patients with G1 grade had a very good prognosis and patients with G2 grade had a relatively poor prognosis, we further divided patients with G2 grade into G2(GIL) and G2(GIH) based on genome instability and found that patients in group G2(GIL) showed significantly longer RFS compared with patients in group G2(GIH) (Figure S2H). Receiver operating characteristic curve (ROC) also showed that this new classification (G1, G2(GIL), G2(GIH), G3) was better than that classification based on genomic instability (GIL, GIH) and histological grading (G1, G2, G3) (Figure S2I). This indicated that the combination of genomic instability and histological grading was of great value in evaluating clinical prognosis. Furthermore, in the new classification system, the proportion of CNVhigh tumor cells in the subcluster was positively correlated with the proportion of tumors from patients with G1 grade and G2(GIH), and was negatively correlated with the proportion of tumors from patients with G2(GIL) (Figure S2J).

We further analyzed the genes specifically expressed in group PHM and found that SMOC1 and PCSK1 were significantly expressed in group PHM (Figure S3A). Indeed, SMOC1 and PCSK1 were reported to be associated with hepatic metastases in NF-Pan-NET.17 Immunohistochemistry assay was performed to detect SMOC1, PCSK1, and PCSK2. In group PHM, the expression of SMOC1 and PCSK1 was significantly increased (Figures S3B and S3C).

Identifying subtype of fibroblasts associated with malignant progression

Fibroblasts was analyzed and visualized in the tSNE map. There were 9 sub-clusters of fibroblasts (Figure 4A). By analyzing the distribution of fibroblasts, we found that sub-cluster 6 was unique to the tumor tissue and was mainly concentrated in group PHM, which was a typical tumor-associated fibroblast (CAFs) (Figure 4B). Sub-cluster 6 particularly expressed NSG1 and COL9A1 (Figure 4C). Sub-cluster 2 was large distributed in group G2(GIH) and sub-cluster 2 particularly expressed LAMP5 and RRGS5 (Figures 4B and 4C). The result of CytoTRACE showed that sub-clusters 2 and 6 had lowest degree of cell differentiation (Figure 4D). Pathway analysis showed that oxidative phosphorylation and pathways of neurodegeneration were mainly enriched in sub-cluster 2 while focal adhesion and PI3K-Akt signaling pathway were mainly enriched in sub-cluster 6 (Figure 4E). Gene-module analysis showed that sub-cluster 6 had activated hippo signaling pathway, cytokine-cytokine receptor interaction, and transforming growth factor (TGF-beta) signaling pathway (Figure 4F). QuSAGE analysis also showed that sub-cluster 6 exhibited activated hedgehog signaling, PI3K-Akt-mTOR signaling pathway, IL6-JAK-STAT3 signaling, IL2-STAT5 signaling, and so on (Figure 4G). By analyzing interleukin family gene expression, we found that sub-cluster 6 might secrete IL32 (Figure 4H). CAF-derived IL32 has been reported to promote cancer cell invasion and metastasis in breast cancer.18 These data indicate that fibroblasts in the tumor tissue were highly heterogeneous and NSG1+ CAFs might play a facilitating role in the metastasis of NF-Pan-NET.Figure 4 Identifying subtype of fibroblasts associated with malignant progression

(A) tSNE map showed the sub-cluster of fibroblasts.

(B) The proportion of sub-cluster in different groups.

(C) Selected marker genes of fibroblasts in cluster 2 and 6.

(D) Differentiated degree of fibroblasts in each cluster was analyzed by CytoTRACE.

(E) Pathway analysis of fibroblasts in cluster 2 and 6.

(F) Gene-module analysis of fibroblasts in each cluster.

(G) The activation of KEGG pathway in fibroblasts in each cluster was analyzed by QuSAGE.

(H) Part interleukin family gene expression of fibroblasts in each cluster.

Analysis of endothelial cell heterogeneity based on genome instability and histological grading

Endothelial cells have been reported to be associated with tumor cell invasion, metastasis, and immunosuppression.19,20 We thus analyzed the subpopulations of endothelial cells and showed in the tSNE map (Figure 5A). According to the expression of marker gene, these 11 clusters of cells were divided into five subtypes (Figure 5B). The proportion of each endothelial cell subpopulation in different groups was calculated and we found that in group PHM the proportion of tip cells and lymphatic endothelial cells were higher. Compared with group G2(GIL), the proportion of tip cells was significantly increased while the proportion of artery cells was significantly decreased in group G2(GIH) (Figure 5C). Gene-module analysis showed that Notch signaling pathway, platelet activation, and PI3K-Akt signaling pathway were enriched in tips cell; Vascular endothelial growth factor (VEGF) signaling pathway and mitogen-activated protein kinase (MAPK) signaling pathway were enriched in artery cells and extracellular matrix (ECM) receptor interaction, cell adhesion molecules and adenosine 5‘-monophosphate (AMP)-activated protein kinase (AMPK) signaling pathway were enriched in lymphatic endothelial cells (Figure 5D). Further analysis showed that CCL21 was significantly expressed in lymphatic endothelial cells (Figure 5E). CCL21 in the tumor microenvironment was reported to be associated with the formation of immune tolerance in cancer.21 To validate this, we performed immunohistochemistry assay to detect PDPN and CCL21 expression in tumor tissues. In group PHM, PDPN+, or CCL21+ endothelial cells were significantly increased (Figures 5F and 5G). In addition, PDPN+ endothelial cells tend to accompanied by CCL21+ endothelial cells (Figure 5H). Both PDPN+ and CCL21+ endothelial cells predicted a worse progress free survival (PFS) probability (Figure 5I). These results suggest that endothelial cells, especially lymphatic endothelial cells, might play key role in the metastasis of NF-Pan-NET.Figure 5 Analysis of endothelial cell heterogeneity based on genome instability and histological grading

(A) tSNE map showed the subpopulations of endothelial cells.

(B) Selected marker genes of each endothelial cell subpopulation.

(C) The proportion of each endothelial cell subpopulation in different groups.

(D) Gene-module analysis of endothelial cells in each cluster.

(E) The expression of chemokines in each endothelial cell subpopulation.

(F) Immunohistochemistry assay showed the typical staining of PDPN and CCL21 in a sample of PHM (primary lesion).

(G) The proportion of PDPN+/CCL21+ endothelial cells in different groups.

(H) The Chi-square test compared the composition ratio of CCL21+ to PDPN+ endothelial cells.

(I) Kaplan-Meier depicted the effect of PDPN+/CCL21+ endothelial cells on PFS probability.

The heterogeneity of immune cells based on genome instability and histological grading

Tumor immune microenvironment plays a vital role in the occurrence and development of tumors and tumor immunotherapy.19,22 We analyzed the subpopulation of monocytic cells and visualized them in the tSNE map (Figures 6A and 6B). We found that macrophages were the most abundant monocytic cells in the microenvironment. CD1A + dendritic cell was significantly enriched in group PHM. CD14+ monocyte was significantly enriched in group G1(GIH) and G2(GIL). CD1C+ dendritic cell as significantly enriched in group G2(GIH) (Figure 6C). We further analyzed the subtypes of macrophage and there were 8 sub-cluster macrophages. Compared with group G2(GIL), the proportion of clusters 2 and 6 were significantly higher and clusters 0 and 4 were significantly lower in group G2(GIH). Cluster 1 were significantly increased group PHM (Figure 6D). Gene-module analysis showed that MAPK signaling pathway, tumor necrosis factor (TNF) signaling pathway, IL-17 signaling pathway were enriched in cluster 1. NF-κB signaling pathway was enriched in Cluster 0. cGMP-PKG signaling pathway and AMPK signaling pathway was enriched in cluster 2 and 6 (Figure 6E). By analyzing interleukin family gene expression, we found that IL-1b, which is an important inducer of immunosuppressive microenvironment in cancer,23 was significantly in clusters 0, 1, and 7 (Figure 6F).Figure 6 The heterogeneity of monocytic cells based on genome instability and histological grading

(A) Selected marker genes of monocytic cells.

(B) tSNE map showed the subpopulations of monocytic cells.

(C) The proportion of subpopulations of monocytic cells in different groups.

(D) The proportion of subpopulations of macrophages in different groups.

(E) Gene-module analysis of macrophages in each cluster.

(F) The expression of IL-1b in each cluster of macrophages.

There were 3,051 CD4+ T cells in our data. The subpopulation of CD4+ T cells were shown in the tSNE map (Figure S4A). The marker genes of CD4+ T cells were was shown in the tSNE map (Figure S4A). There were 4,422 CD8+ T cells in our data. The subpopulation of CD8+ T cells were shown in the tSNE map (Figure S4B). The marker genes of CD8+ T cells were was shown in the tSNE map (Figure S4B). There were 2,889 natural killer cells (NK) and innate lymphoid cells (ILC) in our data. The subpopulation of NK and ILC cells and marker genes were shown in the tSNE map (Figure S4C). We found that Treg cells, mucosal-associated invariant T (MAIT), and ILC were significantly enriched in group PHM (Figure S4D). By multiplex immunohistochemistry, antibody against CD4 was used to stain CD4+ T cells and antibodies against IL2RA and FOXP3 was used to stain Treg cells. We confirmed that the proportion of Treg cells in CD4+ T cells was higher in group PHM (Figures S4E and S4F). Antibody against CD8 was used to stain CD8+ T cells and antibodies against NCR3 and SLC4A10 was used to stain MAIT cells. We found that MAIT was significantly enriched in group PHM (Figures S4G and S4H). These data imply that the PHM had a unique immune cell infiltration pattern.

Discussion

Pan-NENs are a kind of highly heterogeneous tumors. The heterogeneity can be found from clinical symptoms to molecular pathology and distinct prognosis.24 Depending on the degree of differentiation, Pan-NENs can further classified as Pan-NET and Pan-NECs. Pan-NECs are poorly differentiated and highly aggressive. As a comparison, Pan-NET is well-differentiated and indolent.25 According to hormone secretion and the corresponding clinical symptoms, Pan-NET can be divided into functional and non-functional categories.26 Based on the mitotic image and ki-67 proliferation index, Pan-NET are further divided into three grades, G1 (Ki-67 < 2%), G2 (2% ≤ Ki-67 ≤ 20%), and G3(Ki-67 ≥ 20%). This classification forms the basis of current treatment guidelines.4,27 However, the histological grading is dynamic as the disease progresses.8 Furthermore, 90% of Pan-NET are grade G1 or G2, thereby making it hard to decode the malignant potential of Pan-NET.

Attempts to classify the heterogeneity of Pan-NET have never been interrupted. One study involving the WGS data of 98 Pan-NET showed that Pan-NET can be divided into four discrete groups according to CNV patterns. In addition, their data reported an interesting phenomenon that a third of the population possessed limited copy number events and lowest mutation burden.10 This was verified in our present study. We found a certain proportion of endocrine cells could not be identified as malignant by inferCNV due to the low copy number variants. Coincidentally, another study involving the WGS/WES data of 211 Pan-NET gave a similar conclusion: CNV alteration is a better index for NET classification besides histological grading. In their data, there were also approximately one-third of patients without CNV-amplification/deletion. In addition, the distribution of histological grading also varied from the CNV patterns: in the CNV-neutral group, G1 accounted for more than 50%, while in the CNV- amplification/deletion groups, G2 dominated.12

In this study, we performed single-cell sequencing of 17 NF-Pan-NET and found that inferCNV, which is widely used to identify malignant cells in other cancers, had significant deficiencies in identifying tumor cells in scRNA data of NF-Pan-NET. Because a third of the NF-Pan-NET possessed limited copy number variants in bulk-sequences. The proportion of normal endocrine cells in pancreas is very low. Thus, it is reasonable to speculate that endocrine cells in the tumor samples are tumor cells. We further proved that endocrine cells in the tumor samples had tumor cell characteristics by combining scRNA-seq data with that of multi-omics. There was strong heterogeneity between samples. We found that genome instability was better than histological grading in subtyping tumor cells. Particularly, patients with G2 grade can be divided into two separated subgroups based on genome instability. In our single-cell data, clusters with G1 had high CNV and genomic instability. However, in the single-cell data, we only had 4 patients with G1 grade, which may result in an error due to a small sample size. We further perform WES on 77 NF-Pan-NET tumor specimens and analyzed the CNV and genomic instability alteration and found that high or low levels of CNV and genomic instability are present in both G1 and G2. However, in a larger cohort study by Hong et al.,12 the distribution of CNV and genomic instability between G1 and G2 is different. Patients with G1 tend to have lower genomic instability, while patients with G2 tend to have higher genomic instability. Therefore, we focused on the differences of tumor cells under different genome instability patterns, and the differences of tumor immune microenvironment between genome instability and histological grading classification. CAFs affect many aspects of tumor biology, ranging from collagen deposition to immunosuppression, and CAFs have been the target of preclinical and clinical studies due to their potential pro- and anti-tumorigenic functions.28 We found a special subpopulation of fibroblasts that only exist in tumor tissues. These CAFs had an enriched expression of IL32, which was reported to promote cancer cell invasion and metastasis in other cancers implying that CAFs might play a promotional role in the metastasis of NF-Pan-NET. We further analyzed the heterogeneity of endothelial cells under genome instability pattern and histological grading. Compared with group G2(GIL), the proportion of tip cells was significantly increased while the proportion of artery cells was significantly decreased in group G2(GIH). In addition, we found that tip cells and lymphatic endothelial cells might play a significant role in the metastasis of NF-Pan-NET. In terms of the tumor immune microenvironment, we found that there was significant heterogeneity in the enrichment of immune cells regardless of genome instability or histological grading in Pan-NET. Macrophages, CD1A + dendritic cell, Treg, MAIT, and ILC might participate in the process of hepatic metastases in NF-Pan-NET. In summary, we explored the heterogeneity in NF-Pan-NET of different histological grading and genome instability patterns as well as identified the characteristics of PHM. Our results provide valuable insights into better management of NF-Pan-NET.

Conclusion

In brief, we deciphered the heterogeneity of the tumor microenvironment under the pattern of genomic instability and histological grading. In particular, we found that the combination of genomic instability and histological grading was of great value in evaluating clinical prognosis, especially for patients with G2 grade. We also identified the factors associated with hepatic metastases, which might facilitate the understanding of the patterns to decode the malignant potential and of NF-Pan-NET.

Limitation of the study

The number of cases recruited in our study was limited, especially only one patient with liver metastasis. In addition, we had no further experiments to validate our results. We will explore the underlying molecular mechanisms to reveal how genomic instability affects tumor progression in a larger cohort.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Shunrong Ji (jishunrong@fudanpci.org).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• Single-cell RNA-seq data have been deposited at GEO and are publicly available as of the date of publication. Accession numbers are listed in the key resources table.

• This paper does not report original code.

• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Acknowledgments

This work was funded by grants from the 10.13039/501100001809 National Natural Science Foundation of China (81872915 , 82073904 , 82121005 , 81973373 , 21704064 , 82141129 , 82141104 , U21A20374 , and 82303943 ), National Science & Technology Major Project of China—Key New Drug Creation and Manufacturing Program (2018ZX09735–001 , 2018ZX09711002–002–005 ), National Key Basic Research Program of China (2018YFA0507000 ), Shanghai Municipal Education Commission-Shanghai Top-Level University Capacity Building Program (DGF817029-04 ), Hainan Provincial Major Science and Technology Project (ZDKJ2021028 ), Shanghai Municipality Science and Technology Development Fund (21JC1401600 ), Shanghai Municipal Science and Technology Major Project (21JC1401500 ), Scientific Innovation Project of Shanghai Education Committee (2019-01-07-00-07-E00057 ), Clinical Research Plan of Shanghai Hospital Development Center (SHDC2020CR1006A ), Xuhui District Artificial Intelligence Medical Hospital Cooperation Project (2021-011 ). Sailing Project of Science and Technology Commission of Shanghai Municipality (23YF1406900 ). 10.13039/501100002858 China Postdoctoral Science Foundation (2023M730670 ).

Author contributions

Conceptualization: M.-W.W., Y.X., and S.J. methodology: D.Y., X.L., J.G., H.Z., Z.X., and Y.H. investigation: H.Z., Y.Q., J.C., and X.X. visualization: Z.Y. and W.Z. supervision: M.-W.W., Y.X., and S.J. writing—original draft: Z.Y., Q.L., and C.X. writing—review and editing: Z.Y. and C.X. All authors have read and reviewed the manuscript.

Declaration of interests

The authors declare that they have no conflict of interest.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Antibodies	
	
Recombinant Anti-CD4 antibody	Abcam	Cat# ab133616; RRID:AB_2750883	
Recombinant Anti-FOXP3 antibody	Abcam	Cat# ab20034; RRID:AB_445284	
Recombinant Anti-IL-2 Receptor alpha (IL2RA)	Abcam	Cat# ab128955; RRID:AB_11141054	
Recombinant Anti-CD8 alpha antibody	Abcam	Cat# ab237709; RRID:AB_2892677	
SLC4A8/SLC4A10 Polyclonal Antibody	Invitrogen	Cat# PA5-101904; RRID:AB_2851336	
NKp30/NCR3 Polyclonal Antibody	Invitrogen	Cat# PA5-104376; RRID:AB_2853684	
Recombinant Anti-Podoplanin (PDPN) antibody	Abcam	Cat# ab236529; RRID:AB_2889051	
CCL21 Rabbit pAb	Abclonal	Cat# A1896; RRID:AB_2763928	
Recombinant Anti-PC1/3 (PCSK1) Antibody	Abcam	Cat# ab233397; RRID:AB_3096426	
Recombinant Anti-PCSK2 Antibody	Abcam	Cat# ab274418; RRID:AB_3271522	
Recombinant Anti-SMOC1 Antibody	Abcam	Cat# ab313569; RRID:AB_2941846	
Recombinant Anti-SCGN/Secretagogin Antibody	Abcam	Cat# ab137017; RRID:AB_3271547	
	
Biological samples	
	
Pan-NET primary tumor tissue	Fudan University Shanghai Cancer Center	N/A	
Pan-NET adjacent normal tissue	Fudan University Shanghai Cancer Center	N/A	
Pan-NET hepatic metastases	Fudan University Shanghai Cancer Center	N/A	
	
Deposited data	
	
Single RNA-seq of Pan-NET tissue	This paper	GEO: GSE256136	
	
Software and algorithms	
	
GraphPad Prism 8	GraphPad	N/A	
SPSS 26.0	IBM	N/A	

Experimental model and study participant details

Clinical specimens

A total of 17 consecutive patients who underwent surgical treatment at the Fudan University Shanghai Cancer Center (FUSCC) between 2020 and 2021 and were pathologically confirmed to have NF-Pan-NET were recruited for single-cell sequencing in this study. Validation cohort of NF-Pan-NET patients (77) for WES sequencing underwent surgical treatment and were pathologically confirmed to NF-Pan-NET at FUSCC between 2012 and 2021. Patient information and clinicopathological characteristics were acquired from medical records. Progression-free survival was measured from the time after surgery until detection of disease progression or the last follow-up.

Ethics statement

The study was approved by the Ethics Committee of Fudan University Shanghai Cancer Center (2105235-9), and written informed consent was obtained from each patient.

Method details

Single-cell dissociation

During surgery, immediately after the tissue specimens were isolated, the tissue specimens were cuted and kept in MACS Tissue Storage Solution (Miltenyi Biotec) for further processing. Tissue specimens were processed as described below. Briefly, specimen was first washed with phosphate-buffered saline (PBS), minced into small pieces (approximately 1 mm3) on ice and enzymatically digested with collagenase I (Worthington) and DNase I (Worthington) for 45 min at 37°C, with agitation. After digestion, specimen was sieved through a 70 μm cell strainer, and centrifuged at 300 g for 5 min. After the supernatant was removed, the pelleted cells were suspended in red blood cell lysis buffer (Miltenyi Biotec) to lyse red blood cells. After washing with PBS containing 0.04% BSA, the cell pellets were re-suspended in PBS containing 0.04% BSA and re-filtered through a 35μm cell strainer. Dissociated single cells were then stained with AO/PI for viability assessment using Countstar Fluorescence Cell Analyzer. The single-cell suspension was further enriched with a MACS dead cell removal kit (Miltenyi Biotec).

Single-cell sequencing

The scRNA-seq libraries and V(D)J libraries were generated using the 10x Genomics Chromium Controller Instrument and Chromium Single Cell 5’ library & gel bead kit, along with the V(D)J enrichment kit (10x Genomics). Briefly, cells were concentrated to 1000 cells/μL and approximately 10,000 cells were loaded into each channel to generate single-cell Gel Bead-In-Emulsions (GEMs), which resulted in expected mRNA barcoding of 6,000 single-cells for each specimen. After the reverse transcription step, GEMs were broken and barcoded-cDNA was purified and amplified. The amplified barcoded cDNA was used to construct 5’ gene expression libraries and TCR and BCR enriched libraries. For 5’ library construction, the amplified barcoded cDNA was fragmented, A-tailed, ligated with adaptors and index PCR amplified. The final libraries were quantified using the Qubit High Sensitivity DNA assay (Thermo Fisher Scientific) and the size distribution of the libraries were determined using a High Sensitivity DNA chip on a Bioanalyzer 2200 (Agilent). All libraries were sequenced by a NextSeq2000 sequencer (Illumina) on a 150 bp paired-end run.

QuSAGE analysis (gene enrichment analysis)

To characterize the relative activation of a given gene set such as pathway activation, “Human_41GeneSets_190716.gmt” as described before, we performed QuSAGE (2.16.1) analysis.

Differential gene expression analysis

To identify differentially expressed genes among specimen, the function FindMarkers with wilcox rank sum test algorithm was used under following criteria: 1) Log2FC > 0.25; 2) pvalue<0.05; and 3) min.pct>0.1.

CNV estimation

Cells defined as endothelia, fibroblast and SMC were used as reference to identify somatic copy number variations with the R package inferCNV (v0.8.2). We scored each cell for the extent of CNV signal, defined as the mean of squares of CNV values across the genome. CNV correlation was the Pearson correlation between CNV profile of each cell and the average CNV of all malignant cells identified by unsupervised clustering from the same sample.29 Putative CNV-high malignant cells were then defined as those with CNV signal above 0.05 and CNV correlation above 0.5.

Go analysis

Gene ontology (GO) analysis was performed to facilitate elucidating the biological implications of marker genes and differentially expressed genes. We downloaded the GO annotations from NCBI (http://www.ncbi.nlm.nih.gov/), UniProt (http://www.uniprot.org/) and the Gene Ontology (http://www.geneontology.org/). Fisher’s exact test was applied to identify the significant GO categories and false discovery rate (FDR) was used to correct the p-values.

Pathway analysis

Pathway analysis was used to find out the significant pathway of the marker genes and differentially expressed genes according to KEGG database. We turned to the Fisher’s exact test to select the significant pathway, and the threshold of significance was defined by P-value and FDR.

ssGSEA analysis

The ssGSEA function of GSVA package was utilized to quantify the enrichment scores of the geneset for each specimen cells.

Immunohistochemistry (IHC) and multiplex immunohistochemistry (mIHC)

The expression of CD4, FOXP3, IL2RA, CD8, SLC4A10 and NCR3 were detected by multiplex IHC through a tyramide signal amplification system. In brief, slices were deparaffinized with gradient concentrations of xylene and ethanol. Then these slices were incubated with a 3% solution of H2O2 for 15 minutes to suppress endogenous peroxidase activity. Antigen retrieval was performed by heating the samples in the recommended antigen retrieval solution. After blocking with 5% BSA for 30 minutes, the slides were incubated with primary antibodies at optimal dilutions overnight at 4°C and then incubated at 37°C for 1 hours. Next, the slides were incubated with tyramide conjugated fluorophore at room temperature for 20 minutes. The same procedure was repeated with other primary antibodies, and the nuclei were stained with DAPI. The expression of PDPN, CCL21, PCSK1, PCSK2, SMOC1 and SCGN were detected by IHC. The procedure was similar to mIHC. Secondary antibodies were conjugated with enzymes and diaminobenzidine as the chromogen and Mayer’s hematoxylin as the counterstain.

WES, RNA sequencing and proteomic analysis

The experimental procedure was consistent with that of Professor Zhou's previous multi-omics analysis of HBV-related hepatocellular carcinoma.30 Tissue samples were collected within 30 min after operation and snap-frozen in liquid nitrogen. Tissue samples were pulverized using the CryoPrepTM CP02 (Covaris) and then divided into three parts to perform whole exome sequencing (WES), RNA sequencing, proteomic and phosphoproteomic analyses. The specific experimental procedures of each omics were the same as Zhou's research.30

Quantification and statistical analysis

scRNA-seq data analysis was performed by NovelBio Co.,Ltd. with NovelBrain Cloud Analysis Platform (https://www.novelbrain.com). We applied the down specimen analysis among specimen sequenced according to the mapped barcoded reads per cell of each specimen and finally achieved the aggregated matrix. Cells contained over 200 expressed genes and mitochondria UMI rate below 40% passed the cell quality filtering and mitochondria genes were removed in the expression table.

Seurat package (version: 4.0.3, https://satijalab.org/seurat/) was used for cell normalization and regression based on the expression table according to the UMI counts of each specimen and percent of mitochondria rate to obtain the scaled data. Principle component analysis (PCA) was constructed based on the scaled data with top 2 000 high variable genes and top 10 principals were used for t-distributed Stochastic Neighbor Embedding (tSNE) construction and Uniform Manifold Approximation and Projection (UMAP) construction. Utilizing graph-based cluster method, we acquired the unsupervised cell cluster result based the PCA top 10 principal and calculated the marker genes by FindAllMarkers function with wilcox rank sum test algorithm under following criteria: 1) Log2FC> 0.25; 2) P value<0.05; and 3) min.pct >0.1. In order to identify the cell type detailed, the clusters of same cell type were selected for re-tSNE analysis, graph-based clustering and marker analysis. Statistical analysis about the proportion of cell subpopulations was performed by GraphPad Prism. The Student’s t-test was applied for statistical discrepancy between two groups. Kaplan-meier analysis was performed to depict overall survival of patients and the log-rank test was conducted for statistical difference comparisons. P<0.05 was considered significant. ns P>= 0.05, ∗ P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001.

Supplemental information

Document S1. Figures S1–S4

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110836.
==== Refs
References

1 Franko J. Feng W. Yip L. Genovese E. Moser A.J. Non-functional neuroendocrine carcinoma of the pancreas: incidence, tumor biology, and outcomes in 2,158 patients J. Gastrointest. Surg. 14 2010 541 548 10.1007/s11605-009-1115-0 19997980
2 Chauhan A. Kohn E. Del Rivero J. Neuroendocrine Tumors-Less Well Known, Often Misunderstood, and Rapidly Growing in Incidence JAMA Oncol. 6 2020 21 22 10.1001/jamaoncol.2019.4568 31697337
3 Dasari A. Shen C. Halperin D. Zhao B. Zhou S. Xu Y. Shih T. Yao J.C. Trends in the Incidence, Prevalence, and Survival Outcomes in Patients With Neuroendocrine Tumors in the United States JAMA Oncol. 3 2017 1335 1342 10.1001/jamaoncol.2017.0589 28448665
4 Pavel M. Öberg K. Falconi M. Krenning E.P. Sundin A. Perren A. Berruti A. ESMO Guidelines Committee Electronic address clinicalguidelines@esmoorg Gastroenteropancreatic neuroendocrine neoplasms: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up Ann. Oncol. 31 2020 844 860 10.1016/j.annonc.2020.03.304 32272208
5 Frilling A. Modlin I.M. Kidd M. Russell C. Breitenstein S. Salem R. Kwekkeboom D. Lau W.Y. Klersy C. Vilgrain V. Recommendations for management of patients with neuroendocrine liver metastases Lancet Oncol. 15 2014 e8 e21 10.1016/S1470-2045(13)70362-0 24384494
6 Raymond E. Dahan L. Raoul J.L. Bang Y.J. Borbath I. Lombard-Bohas C. Valle J. Metrakos P. Smith D. Vinik A. Sunitinib malate for the treatment of pancreatic neuroendocrine tumors N. Engl. J. Med. 364 2011 501 513 10.1056/NEJMoa1003825 21306237
7 Yao J.C. Hassan M. Phan A. Dagohoy C. Leary C. Mares J.E. Abdalla E.K. Fleming J.B. Vauthey J.N. Rashid A. Evans D.B. One hundred years after "carcinoid": epidemiology of and prognostic factors for neuroendocrine tumors in 35,825 cases in the United States J. Clin. Oncol. 26 2008 3063 3072 10.1200/JCO.2007.15.4377 18565894
8 Grillo F. Albertelli M. Brisigotti M.P. Borra T. Boschetti M. Fiocca R. Ferone D. Mastracci L. Grade Increases in Gastroenteropancreatic Neuroendocrine Tumor Metastases Compared to the Primary Tumor Neuroendocrinology 103 2016 452 459 10.1159/000439434 26337010
9 Keck K.J. Choi A. Maxwell J.E. Li G. O'Dorisio T.M. Breheny P. Bellizzi A.M. Howe J.R. Increased Grade in Neuroendocrine Tumor Metastases Negatively Impacts Survival Ann. Surg Oncol. 24 2017 2206 2212 10.1245/s10434-017-5899-y 28560597
10 Scarpa A. Chang D.K. Nones K. Corbo V. Patch A.M. Bailey P. Lawlor R.T. Johns A.L. Miller D.K. Mafficini A. Whole-genome landscape of pancreatic neuroendocrine tumours Nature 543 2017 65 71 10.1038/nature21063 28199314
11 Ciriello G. Miller M.L. Aksoy B.A. Senbabaoglu Y. Schultz N. Sander C. Emerging landscape of oncogenic signatures across human cancers Nat. Genet. 45 2013 1127 1133 10.1038/ng.2762 24071851
12 Hong X. Qiao S. Li F. Wang W. Jiang R. Wu H. Chen H. Liu L. Peng J. Wang J. Whole-genome sequencing reveals distinct genetic bases for insulinomas and non-functional pancreatic neuroendocrine tumours: leading to a new classification system Gut 69 2020 877 887 10.1136/gutjnl-2018-317233 31462556
13 Ben-David U. Amon A. Context is everything: aneuploidy in cancer Nat. Rev. Genet. 21 2020 44 62 10.1038/s41576-019-0171-x 31548659
14 Turajlic S. Sottoriva A. Graham T. Swanton C. Resolving genetic heterogeneity in cancer Nat. Rev. Genet. 20 2019 404 416 10.1038/s41576-019-0114-6 30918367
15 Peng J. Sun B.F. Chen C.Y. Zhou J.Y. Chen Y.S. Chen H. Liu L. Huang D. Jiang J. Cui G.S. Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in pancreatic ductal adenocarcinoma Cell Res. 29 2019 725 738 10.1038/s41422-019-0195-y 31273297
16 Alvarez M.J. Subramaniam P.S. Tang L.H. Grunn A. Aburi M. Rieckhof G. Komissarova E.V. Hagan E.A. Bodei L. Clemons P.A. A precision oncology approach to the pharmacological targeting of mechanistic dependencies in neuroendocrine tumors Nat. Genet. 50 2018 979 989 10.1038/s41588-018-0138-4 29915428
17 Zhou Y. Liu S. Liu C. Yang J. Lin Q. Zheng S. Chen C. Zhou Q. Chen R. Single-cell RNA sequencing reveals spatiotemporal heterogeneity and malignant progression in pancreatic neuroendocrine tumor Int. J. Biol. Sci. 17 2021 3760 3775 10.7150/ijbs.61717 34671197
18 Wen S. Hou Y. Fu L. Xi L. Yang D. Zhao M. Qin Y. Sun K. Teng Y. Liu M. Cancer-associated fibroblast (CAF)-derived IL32 promotes breast cancer cell invasion and metastasis via integrin beta3-p38 MAPK signalling Cancer Lett. 442 2019 320 332 10.1016/j.canlet.2018.10.015 30391782
19 de Visser K.E. Joyce J.A. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth Cancer Cell 41 2023 374 403 10.1016/j.ccell.2023.02.016 36917948
20 Geldhof V. de Rooij L.P.M.H. Sokol L. Amersfoort J. De Schepper M. Rohlenova K. Hoste G. Vanderstichele A. Delsupehe A.M. Isnaldi E. Single cell atlas identifies lipid-processing and immunomodulatory endothelial cells in healthy and malignant breast Nat. Commun. 13 2022 5511 10.1038/s41467-022-33052-y 36127427
21 Shields J.D. Kourtis I.C. Tomei A.A. Roberts J.M. Swartz M.A. Induction of lymphoidlike stroma and immune escape by tumors that express the chemokine CCL21 Science 328 2010 749 752 10.1126/science.1185837 20339029
22 Chao Y. Liu Z. Biomaterials tools to modulate the tumour microenvironment in immunotherapy Nat. Rev. Bioeng. 1 2023 125 138 10.1038/s44222-022-00004-6
23 Kaplanov I. Carmi Y. Kornetsky R. Shemesh A. Shurin G.V. Shurin M.R. Dinarello C.A. Voronov E. Apte R.N. Blocking IL-1beta reverses the immunosuppression in mouse breast cancer and synergizes with anti-PD-1 for tumor abrogation Proc. Natl. Acad. Sci. USA 116 2019 1361 1369 10.1073/pnas.1812266115 30545915
24 Kidd M. Modlin I. Öberg K. Towards a new classification of gastroenteropancreatic neuroendocrine neoplasms Nat. Rev. Clin. Oncol. 13 2016 691 705 10.1038/nrclinonc.2016.85 27273044
25 Kawasaki K. Rekhtman N. Quintanal-Villalonga Á. Rudin C.M. Neuroendocrine neoplasms of the lung and gastrointestinal system: convergent biology and a path to better therapies Nat. Rev. Clin. Oncol. 20 2023 16 32 10.1038/s41571-022-00696-0 36307533
26 Cives M. Strosberg J.R. Gastroenteropancreatic Neuroendocrine Tumors CA. Cancer J. Clin. 68 2018 471 487 10.3322/caac.21493 30295930
27 Garcia-Carbonero R. Anton-Pascual B. Modrego A. Del Carmen Riesco-Martinez M. Lens-Pardo A. Carretero-Puche C. Rubio-Cuesta B. Soldevilla B. Advances in the Treatment of Gastroenteropancreatic Neuroendocrine Carcinomas: are we moving forward? Endocr. Rev. 44 2023 724 736 10.1210/endrev/bnad006 36879384
28 Caligiuri G. Tuveson D.A. Activated fibroblasts in cancer: Perspectives and challenges Cancer Cell 41 2023 434 449 10.1016/j.ccell.2023.02.015 36917949
29 Zhang L. Li Z. Skrzypczynska K.M. Fang Q. Zhang W. O'Brien S.A. He Y. Wang L. Zhang Q. Kim A. Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer Cell 181 2020 442 459.e29 10.1016/j.cell.2020.03.048 32302573
30 Gao Q. Zhu H. Dong L. Shi W. Chen R. Song Z. Huang C. Li J. Dong X. Zhou Y. Integrated Proteogenomic Characterization of HBV-Related Hepatocellular Carcinoma Cell 179 2019 561 577.e22 10.1016/j.cell.2019.08.052 31585088
