Identification of the Immune Subtypes for the Prediction of Metastasis in Pancreatic Neuroendocrine Tumors

Pancreatic neuroendocrine tumors (pNETs) are a subgroup of neuroendocrine neoplasms that have relatively distinct biological behavior and clinical management compared with pancreatic adenocarcinoma [1]. The incidence of pNETs has risen rapidly in recent years, and it has become a significant health burden worldwide [2]. Although most pNETs are inert tumors, they have a high rate of metastasis [3]. Surgical resection is the only curative treatment for localized tumors [4]. Unfortunately, early diagnosis is rare due to a lack of specific clinical symptoms. Most patients are diagnosed in an advanced stage when most systemic therapies lack an objective response, and the disease causes severe morbidity [5]. Therefore, prediction of the risk of metastasis in patients with pNETs is an important clinical concern.

Tumor microenvironment (TME) denotes the noncancerous cells and components presented in the tumor, including molecules produced and released by them [6]. Noncancerous cells in tumors, such as stromal cells and immune cells, have been known for quite some time to play critical roles in tumor initiation, growth, invasion, metastasis, and response to therapies [7]. For instance, Chen et al. [8] found that cancer-associated fibroblast-induced M2-polarized macrophages can promote hepatocellular carcinoma progression via the plasminogen activator inhibitor-1 pathway. In addition, liver metastases siphon activated CD8⁺ T cells from systemic circulation which then undergo apoptosis following their interaction with monocyte-derived macrophages [9]. Furthermore, a recent study demonstrates that granulin secretion by metastasis-associated macrophages activates resident hepatic stellate cells into myofibroblasts that secrete periostin, resulting in a fibrotic microenvironment that sustains metastatic tumor growth [10]. The infiltrating immune cells also play critical roles in pNETs. Greenberg et al. [11] found that the primary tumors of metastatic pNETs showed significant activation of immune-related pathways and metastatic tumors featured increased numbers of tumor-infiltrating T cells compared to localized tumors. A cohort revealed that high mast cell density was an independent predictor of prolonged progression-free survival in pNETs [12]. Besides, Tanno et al. [13] found that pNETs are immune cold with no evidence of T-cell exclusion and the low density of immune infiltrates indicates poor antitumor immune responses using tissue microarrays from 61 pNETs patients. However, compared with other malignancies, the immune infiltration research in pNETs remains sparse.

The arrival of modern immunotherapy has changed the landscape of cancer treatment; prime examples include immune checkpoint inhibitors and chimeric antigen receptor T-cell therapies [14, 16]. Accumulating evidence indicates that immunotherapy seems to be a promising treatment option for patients with pNETs [17]. According to a prospective multicenter clinical trial CA209-538, combination immunotherapy with ipilimumab and nivolumab demonstrated significant clinical activity in patients with advanced pNETs [18]. In the KEYNOTE-028 trial, pembrolizumab demonstrated some antitumor activity in patients with pNETs [19]. However, there is still no consensus on the effect of immunotherapy in pNETs, and studies that explore the relation between immune infiltration and clinical characteristics in patients with pNETs have been largely inadequate.

In this study, we analyzed the infiltration pattern of immune cells within the primary and metastatic tumors of pNETs using the single-sample Gene Set Enrichment Analysis (ssGSEA) method. Three subtypes of different immune infiltration degrees were identified, and we explored their relations with tumor metastasis. Moreover, differentially expressed genes (DEGs) among three clusters were identified and then subjected to enrichment analysis. We believe that our study could assist in predicting the metastasis of pNETs patients and can further deepen our understanding of the TME of pNETs.

We obtained the gene expression profiling dataset GSE73338 generated by human oligo microarrays from the Gene Expression Omnibus (GEO) database. Literature mining was carried to obtain the clinical characteristics of the core clinical samples used in GSE73338, including 75 primary tumor samples, 4 adjacent normal samples, 4 liver metastasis samples, and 3 lymph node metastasis samples [20]. The clinical and pathological characteristics of these patients are summarized in Table 1. We obtained the RNA-Seq data GSE178398 from the GEO database, which contains 22 primary pNET samples (15 with localized tumors and 7 with metastasis).

The ssGSEA algorithm was based on 29 immune gene sets, including genes related to different immune cell types, functions, pathways, and checkpoints [21]. We employed the ssGSEA algorithm via GSVA packages of R to comprehensively assess the immune cell infiltration degrees of tumor samples in the study. The Estimation of Stromal and Immune Cells in Malignant Tumors using the Expression Data (ESTIMATE) algorithm was utilized to calculate the tumor purity, stromal score, immune score, and ESTIMATE score of tumor samples via estimate package of R [22]. Estimating the Proportion of Immune and Cancer cells (EPIC) online software (http://epic.gfellerlab.org/) was utilized to assess the proportion of different immune cells of tumor samples in the study [23]. We compared the results obtained by the ssGSEA and EPIC using Pearson’s correlation. Statistically significant thresholds were set at a p < 0.05.

By applying the unsupervised k-means clustering algorithm, pNET samples were assigned into three clusters based on the ssGSEA scores of infiltrated immune cells. pNET samples in GSE73338 were classified as high immune cell infiltration cluster, middle immune cell infiltration cluster, and low immune cell infiltration cluster using stats package.

The DEGs among three immune cell infiltration groups were identified by the limma package using the unpaired algorithm. The differential expressed genes between primary tumor samples and paired metastatic tumor samples were identified by the limma package using paired algorithm.

The protein-protein interaction network of DEGs among three immune cell infiltration subtypes was constructed using the STRING online database (https://string-db.org/). An interaction score of >0.4 and p < 0.05 were regarded as statistically significant. To understand the underlying biological mechanisms, Gene Ontology (GO) annotation and pathway analyses of DEGs were conducted using STRING, Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/), and Reactome (https://reactome.org/). Statistically significant thresholds were set at a p < 0.05.

All statistical analyses were conducted using the R version 3.6.3 (https://www.r-project.org/). Continuous variables were analyzed using Student’s t test. Categorical variables were analyzed using the χ² test. The relation between immune cell infiltration degrees and clinical characteristics was evaluated using logistic regression analysis. The results with p < 0.05 were regarded statistically significant.

We obtained 75 pNET samples’ data from the GEO database. To further evaluate the subpopulations of inflammatory cells and enrichment of immune signature in pNETs, the ssGSEA algorithms were performed. Then, we obtained the richness levels of 29 immune-related cells in pNET samples grouped by grading and metastasis (shown in Fig. 1). We found that all types of immune cells except immature dendritic cells (DCs) had higher levels of infiltration in the metastatic group. To further confirm our finding, we evaluated the relation between immune cell infiltration degrees and metastasis using logistic regression analysis (shown in Table 2). We found that metastasis and all types of immune cells except immature DCs were positively correlated, which was consistent with our previous finding.

Noticing the relation between immune cell infiltration and metastasis, we further carried out immune subtype analysis. By applying the unsupervised k-means clustering algorithm, pNET samples were assigned into three clusters based on immune infiltration degrees, including cluster 1 (18 samples), cluster 2 (51 samples), and cluster 3 (6 samples). These clusters were termed as low immune cell infiltration subtype (Immunity-L), middle immune cell infiltration subtype (Immunity-M), high immune cell infiltration subtype (Immunity-H), respectively (shown in Fig. 2a, b). To validate the above clustering result, tumor purity, immune score, stromal score, and ESTIMATE score were calculated using ESTIMATE. We found that high immune cell infiltration degrees were positively correlated with the stromal score, immune score, and ESTIMATE score, but tumor purity was opposite (shown in Fig. 2c–f). Besides, immune cell infiltration subtypes could predict metastasis, which was consistent with our previous findings (shown in Fig. 2h). EPIC was utilized to validate the immune cell proportions calculated by ssGSEA. Except for CD8+ T cells, the results obtained by the two methods are consistent (shown in Fig. 3). This inconsistency may due to a different algorithm and classification of immune cells. To further validate the relation between immune cell infiltration and metastasis, a validation cohort (GSE178398) was analyzed in this study, and we obtained consistent results (shown in Fig. 4).

The limma package of R was used to screen for DEGs among three immune subtypes. As shown in Figure 5a, a total of 80 DEGs were obtained. The heat map of the DEGs is shown in Figure 5b. The protein-protein interaction network of DEGs was constructed using STRING. An interaction score of >0.4 and p < 0.05 indicated statistical significance (shown in Fig. 6). Twenty-one of the 80 DEGs were used in ssGSEA and ESTIMATE, including IGSF6, TFEC, TREM2, FCER1G, SLC15A3, CTSS, CCR5, CSTA, CD52, CD37, IFI16, EMP3, ITGAL, LGALS1, SLA, CD48, ITGAM, MNDA, PDCD1LG2, GPSM3, CCL5, DAB2, MS4A6A, CXCR4, IL18BP, RUNX2, GZMK, FCGR1A, and CCR7. Considering their functional correlation with the remaining 51 DEGs, we did not delete these genes in the next functional enrichment analysis. The results of the GO enrichment analysis for the DEGs are shown in Figure 7. In the biological process group, these genes were mainly enriched in immune-related process, such as immune response, leukocyte activation, inflammatory response, and immune effector process. In the molecular function group, genes were mainly enriched in signaling receptor activity, C-C chemokine receptor activity, and cytokine binding. Then, pathway enrichment was performed using KEGG and Reactome. We found that these genes predominantly participated in immune-related pathways, including immune system, neutrophil degranulation, cytokine signaling, and signaling by interleukins.

Noticing the correlation between immune infiltration and metastasis, we further explored the expression of metastasis-related genes in different immune subtypes. A list of 27 metastasis-related genes was obtained from cBioPortal. We found that the expression pattern of these genes was different in three immune subtypes (shown in Fig. 8). Eleven of 27 metastasis-related genes were differentially expressed (Table 3).

We next compared the immune infiltration patterns of three liver and three lymph node metastatic tumor samples with their paired primary tumor samples (shown in Fig. 9). Despite differences in immune infiltration patterns between primary and metastatic tumors, paired samples can still be clustered together except for one pair of samples (primary tumor 2 and hepatic metastasis 2). In addition, we found that liver metastatic tumors had higher tumor purity and a lower ESTIMATE score compared with paired primary pNETs tumors. The immune score was not significantly different in both liver and lymph node metastatic tumors compared with primary tumors. The results are presented in Figure 10. We next carried out a differential analysis between liver metastasis samples and lymph node metastasis samples. However, no significantly differentially expressed gene was identified, and the top 10 genes are shown in Table 4.

Although pNETs are commonly considered to be more indolent than their exocrine counterpart pancreatic ductal adenocarcinoma, they do have the propensity to metastasize to regional lymph nodes and distant metastatic sites [24]. From the clinical point of view, treatment management for metastatic pNETs is very challenging. Although there are still some therapeutic options for patients with metastatic pNETs, response rates and clinical benefits have been modest [25]. An improved understanding of the immune-mediated regulatory mechanisms underlying the metastatic process may provide novel strategies to improve treatment outcomes.

In the present study, we noticed that all types of immune cells except immature DC had higher levels of infiltration in the metastatic group. DCs are professional and the most powerful antigen-presenting cells. Conversion of a DC from an immature to mature DC is important for initiation of immune responses [26]. Compared with the mature DC, the immature DC has reduced expression of activation markers, reduced ability to stimulate T cells, and reduced ability to migrate [27]. These immature cells are key orchestrators of immunosuppressive microenvironment in many malignancies [28]. We found that both immature DC and immunosuppression were negatively correlated with metastasis in pNETs. Besides, the infiltration level of other immune cells such as CD8⁺ T cells and natural killer cells was elevated in primary pNETs tumors with metastasis. Therefore, we can infer that immature DCs can reduce immune cell infiltration by creating an immunosuppressive environment, thereby restraining tumor metastasis, heralding a possible treatment method for patients with pNETs. In the present study, we found more cytotoxic lymphocytes such as CD8⁺ T cells and natural killer cells to be infiltrated in the tumor tissues, and the immune response was higher in primary pNETs tumors with metastasis. However, regulatory T cell which acts as a negative immunomodulator was also found elevated in cases with metastasis. This may seem contradictory. But the immune system is an interactive and complex network involving various kinds of cells, humoral factors, and cytokines. It is hard to evaluate the immune status just by one type of immune cell. In addition, some immune cells have bidirectional effects on tumorigenesis and progression. For instance, IL-17A secreted by Th17 could promote the proliferation of malignant cells and induce angiogenic constituents by stimulating fibroblasts, resulting in tumor progression [29]. It has been reported that Th17 correlated with reduced disease-free survival rates in many malignancies, such as pancreatic cancer and breast cancer [30, 31]. On the other hand, Th17 showed tumor-suppressing effect by promoting inflammation [32]. Further investigations are needed to elucidate the exact role of these cells in pNETs. Then, we identified three subtypes with different immune cell infiltration degrees by applying the unsupervised k-means clustering algorithm. We found that immune cell infiltration subtypes could predict metastasis. At present, only a few articles about immune infiltration in pNETs can be found. Zhang et al. [33] explored the infiltration degrees of 9 immune cell types by immunohistochemistry and discovered that tumor-infiltrating neutrophils were independent factors for poor overall survival and relapse-free survival. Debien et al. [34] identified two pNET subgroups based on the neutrophils signature and complement pathway genes and revealed the involvement of neutrophils in metastatic evolution of pNETs. Tumor-associated macrophages in pNETs might correlate with tumor progression and metastasis formation [35]. In contrast to previous studies, we analyzed the whole immune infiltration patterns in pNETs using ssGSEA, which can help to obtain a more comprehensive result.

Noncancerous cells in tumors, including stromal cells and immune cells, play critical roles in tumor initiation, growth, invasion, metastasis, and response to therapies [7]. In our study, we obtained higher infiltration scores of immune cells and stromal cells in pNETs with metastasis. In other words, the primary sites of pNETs with metastasis contained more immune cells and stromal cells and less tumor cells. Therefore, tumor purity was negatively correlated with metastasis.

In our study, we found that the expression pattern of metastasis-related genes was different in three immune subtypes. Eleven metastasis-related genes were differentially expressed, including MMP14, MMP2, MMP12, MMP7, SPARC, MMP19, ITGAV, MMP23B, MMP1, MMP25, and MMP9. Matrix metalloproteinases (MMPs) are members of a zinc-dependent endopeptidase family, which plays vital roles in initiation, proliferation, and metastasis of cancer through the breakdown of extracellular matrix physical barriers [36, 39]. Overexpression of MMPs is associated with poor prognosis of cancer [40]. There is increasing evidence implicating MMPs as major regulators of immunity and TME. For instance, CCL5 of glioma-associated microglia/macrophages modulated the migratory and invasive activities of human glioma cells in association with MMP2 expression [41]. Mizuno et al. [42] found that tumor-associated neutrophils contribute to the tumor invasion and angiogenesis through the production of MMP9. However, to the best of our knowledge, no research on immune cell infiltration and MMP has been reported in pNETs. Therefore, further investigations are needed to clarify the immune-related functions of these genes in pNETs. The ESTIMATE algorithm can predict the level of infiltrating stromal and immune cells, and these form the basis for the ESTIMATE score to infer tumor purity in tumor tissue [22]. The ESTIMATE score is calculated by combining the stromal score and immune score, and tumor purity is then predicted based on the ESTIMATE score using the nonlinear least squares method [22]. Although the difference is not significant, we found that the immune score and stroma score were lower in metastatic tumors. But the immune score plus stroma score, which equals to the ESTIMATE score, and the further calculated tumor purity were significantly different between primary and metastatic tumors. Our study indicates that metastatic tumors contain a higher proportion of tumor cells and a lower proportion of noncancerous cells. These noncancerous cells can facilitate desmoplasia and immunosuppression in primary and metastatic sites or can promote metastasis by stimulating angiogenesis, epithelial-mesenchymal transition, invasion, and pre-metastatic niche formation [43]. So far, relevant studies are sparse, and we do not know the cause of the difference in the composition of primary and metastatic tumors and the effect of this difference on disease progression. Large-scale multicenter clinical samples and further fundamental studies are required.

Several limitations to the present study should be mentioned. First, lack of large-scale multicenter clinical samples and some other clinical pathological features, such as treatment and survival information. Besides, we only know whether metastasis was present at the time of diagnosis. Observation period data are absent in our research. Therefore, we did not explore the correlation between the immune infiltration with immunotherapy responsiveness and prognosis. Second, we found that immune cell infiltration subtypes could predict metastasis, and the expression pattern of metastasis-related genes was different in these subtypes. However, the research of underlying regulatory mechanism was not included and need to be studied further. Third, a lack of biological validation is present in our study.

In conclusion, our study systematically evaluated the landscape of immune cells in pNET samples and identified three immune cell infiltration subtypes. We found that these subtypes could predict metastasis, and functional enrichment analysis of characteristic genes was carried out. Then, we identified 11 differentially expressed metastasis-related genes among these subtypes and compared the immune infiltration patterns between primary and metastatic tumors. Our findings may deepen the understanding of the immune-mediated regulatory mechanisms underlying pNETs and may provide some promising targets for immunotherapy.

We thank BioMed Proofreading (http://www.biomedproofreading.com/) for its linguistic assistance during the preparation of this manuscript.

Ethical approval and consent were not required as this study was based on publicly available data.

The authors declare no conflicts of interest for the publication of this manuscript.

The study was supported by the Educational Commission of Liaoning Province, People’s Republic of China (Grant No. L2014294).

The authors contributed to this study and manuscript in the following manner: data collection – D.M. and L.Z.; statistical analysis – C.Z. and D.M.; writing and editing – C.Z. and J.L.; supervision – C.G.; funding acquisition – C.G. All the authors read and approved the final manuscript.

The data used to support the findings of this study are included within the article. Further inquiries can be directed to the corresponding author.