Introduction

Neuroendocrine cells are present throughout the body and are thought to give rise to neuroendocrine tumors (NETs) in multiple organs. Small intestine NETs (SiNETs) have the second highest incidence, after lung NETs [1]. Histopathological grading of SiNET is based on mitoses rate and on the Ki67 proliferation index, with low-grade comprising G1 and G2, and high-grade comprising G3. In addition, neuroendocrine carcinomas (NECs) are included in the G3 category, and have poorly differentiated histopathological appearance [2].

The majority of SiNETs are low grade well-differentiated tumors, usually diagnosed at advanced stages and often present with metastases and stage IV tumor, highlighting their high metastatic potential [3]. Interestingly, SiNETs are often multifocal, with a common finding of synchronous multiple primaries along the small intestine [4, 5]. Two recent studies demonstrated that synchronous SiNETs primaries within the same individual usually do not share somatic alterations, indicating that they arise from different cell clones [6, 7].

SiNETs are usually sporadic, with a surprisingly low mutational burden or recurrent mutations, and show mainly chromosomal alterations [8-10]. The most frequent chromosomal aberration found in SiNETs is the loss of chromosome 18 as seen in 61-89% of patients [11]. In a cohort of 180 SiNETs, the most frequent mutation (found in 8% of tumors) were heterozygous frameshift mutations of CDKN1B, encoding the p27 cyclin regulator, suggesting a role of cell cycle deregulation in the progression of SiNETs [12]. Further studies reported heterogeneity of the CDKN1B mutations even among lesions of the same patient [13]. Overall, there is a limited frequency and consistency of genetic aberrations in SiNETs, of which mechanisms of initiation and progression remain poorly understood.

Another intriguing phenomenon with respect to the large group of neuroendocrine neoplasms (NENs) is that they are sometimes accompanied by a carcinoma component, resulting in “mixed” neuroendocrine non-neuroendocrine neoplasms (MiNEN) [14]. For a NEN to be termed mixed, it is necessary for each component to comprise at least 30% of the tumor, suggesting that many additional tumors may have both components but are not recognized as MiNEN. The origin of MiNENs remains unknown, and their complexity hinders conventional treatment approaches.

In this study, we sought to investigate the molecular characteristics of SiNETs to uncover novel insights into their oncogenic transformation, the resulting tumor ecosystems, and potential vulnerabilities. To this end, we applied single-cell transcriptomics to SiNETs a well as to one MiNEN.

Results

Single cell and single nuclei RNA-seq profiling of SiNETs

To understand the tumor ecosystem of low-grade SiNETs, we generated scRNA-seq profiles for ten primary tumor samples from eight patients using 10x Chromium (Table S1). We profiled single cells from three fresh surgical samples (SiNET1-3) and single nuclei from the remaining seven frozen samples (SiNET4-10). After initial quality controls, we retained 29,198 cells from the ten patients (see Methods). For comparison, we also profiled one fresh small intestine adenocarcinoma (SiAdeno) and included it in our analyses.

We identified clusters of cells in each tumor and annotated them based on the expression of known marker genes as neuroendocrine (NE) cells, and as tumor microenvironment cells (TME) cells including T cells, B/plasma cells, macrophages, fibroblasts, endothelial cells, epithelial cells and natural killer (NK) cells (Table S2). Clustering and marker expression of two exemplary SiNETs are shown in Fig. 1A-D, and the remaining tumors are shown in Fig. S1. Cell type proportions varied considerably between tumors (Fig. 1E). In five tumors NE cells were the most frequent, while other tumors had the highest frequency of other cell types, such as endothelial, fibroblast, epithelial or T cells. The unique composition of each tumor sample likely reflects a combination of tumor-specific biology with spatial sampling within the tumor. Technical effects (e.g. fresh vs. frozen samples) could also impact the capture of distinct cell types, although we did not observe a clear pattern of such bias.

Figure 1.

Heterogeneity of NE cells reveals two SiNET subtypes

We next focused on the NE cells and defined both their common and their variable expression profiles across patients. To this end, we examined the eight samples (one per patient) in which we identified NE cells. We first compared NE cells from each tumor to non-NE cells from the same tumor and identified all genes that are upregulated in NE cells. Genes that were consistently upregulated in NE cells from most tumors (n>5) were defined as common NE genes (Fig. 2A). These signature genes were specifically upregulated in NETs of the small intestine compared to other NETs that were profiled by bulk RNA-seq [15], therefore defining an SiNET signature (Fig. 2B, Table S3). The SiNET signature consists of 26 genes, including CHGA, CHGB, and TPH1, known markers of SiNET [16] as well as of enterochromaffin cells [17]. These genes were most significantly enriched with gene-sets associated with neuroendocrine-related functions such as exocytosis, Neurosecretory vesicle, and serotonergic synapse (all with adjusted P < 0.001).

Figure 2.

However, many more genes were found as upregulated in NE cells of only few of the SiNETs, and these were separated into tumor specific (n=1, see Fig. S2A) and subset-specific (between 2 and 5 tumors). Subset-specific genes may help to uncover functionally distinct subtypes of SiNETs and therefore we investigated them further. These genes exhibited nine distinct expression patterns across the ten SiNETs (Fig. 2C, Table S3). Notably, pattern #1, which was associated with the largest number of genes (n=26), was highly expressed in three tumors that had low expression of all other eight patterns. In contrast, the other seven tumors highly expressed genes from multiple patterns, and had low expression of pattern #1 genes. This analysis highlighted a division of the ten SiNETs into two primary clusters based on pattern #1 vs. all other patterns. Notably, the same division was clearly detected by clustering of the SiNETs based on the global expression profiles of NE cells in each tumor (Fig. 2D).

The first subtype (SiNETs 1, 2 and 9) was associated with upregulation of 129 genes, including the epithelial markers EPCAM and KRT8 (Table S4). The second subtype (SiNET3-8, 10) was associated with upregulation of 73 genes, including neuronal-related genes such as NAV3, RYR2, CACNA1A/B, KCNIP1. This led us to speculate that even though both subtypes express the neuroendocrine markers described above (i.e. common genes), the first subtype may have higher overall similarity to epithelial cells compared to the second subtype, which might have more neuronal features.

To examine this possibility, we evaluated the expression of other cell type signatures across NE cells from the two SiNET subtypes (Fig. 2E). Signatures of multiple epithelial cell types from the small intestine or from the lung, including those of enteroendocrine cells, were expressed more highly in subtype 1 than in subtype 2 SiNETs. Nevertheless, the expression of such epithelial signatures was considerably lower in NE cells from subtype 1 tumors than in bone-fide epithelial cells such as the non-malignant epithelial cells identified in SiNET6 or the malignant cells in the siAdeno tumor. Thus, epithelial signatures are upregulated in subtype 1 tumors but not at the same level as in bone-fide epithelial cells. Conversely, signatures of neuronal cell types had the opposite expression pattern – highest in subtype 2 SiNETs, intermediate in subtype 1 SiNETs and lowest in bone-fide epithelial cells. Based on these results we renamed subtypes 1 and 2 as the epithelial-like and neuronal-like SiNET subtypes.

To validate the existence of these two SiNET subtypes, we turned to analyze an external bulk RNA-seq dataset of SiNETs [15]. The signal of NE cells is diluted in bulk RNA-seq due to the profiling of entire tumor samples, and especially given the wide diversity of cell type compositions that we observe in single cell analysis (see Fig. 1E). However, we wondered whether we would still be able to detect the two subtypes when directly analyzing the subset-specific genes defined above. Indeed, the 81 SiNET bulk profiles could be separated into the epithelial-like and neuronal-like subtypes based on expression of pattern #1 genes vs. the genes of the other patterns (Fig. 2F). Moreover, the proportion of the two subtypes (20% vs. 80%) were comparable to those seen in our single cell cohort (30% vs. 70%).

Heterogeneity in the SiNET tumor microenvironment

Next, we examined the diversity of cellular states within each non-NE cell type of the tumor microenvironment (TME). For each cell type, we analyzed the diversity within each tumor, searching for distinguishable subpopulations of cells, their occurrence across tumors and their potential functional implications (see Fig. 3). Below we briefly describe the three most notable cases of diversity within cell-type that we observed.

Figure 3.

Among fibroblasts, we found the highest diversity in SiNET8, with three distinguishable subpopulations (Fig. 3A). To understand their differences, we both conducted a differential expression analysis (Fig. 3B) and compared their profiles to recently described fibroblast signatures [18] (Fig. 3C). This analysis indicated that cluster 1 resembles other cancer-associated fibroblasts (CAFs), while cluster 2 and 3 are distinguished by specific expression profiles. Cluster 3 cells were primarily distinguished by expression of MHC-II genes and are thus consistent with previous observations of antigen-presenting CAFs [19]. Cluster 2 cells upregulate several programs of pericytes and myofibroblasts, but are further distinguished by upregulation of four ABC-transporters (ABCA6, ABCA8, ABCA9 and ABCA10). To our knowledge, such consistent upregulation of multiple ABC transporters was not observed in previous analysis of CAF heterogeneity [20, 21]. Yet, in our data, upregulation of ABC transporters was also detected in a subset of CAFs from SiNET3 and SiNET5 (Fig. S3A-D). These results suggest a unique feature of CAFs in a subset of SiNETs, possibly due to the unique microenvironment of SiNETs.

Among endothelial cells, we found two highly distinct subpopulations in SiNET5 (Fig. 3D-E). Comparison to pan-cancer endothelial signatures mapped the smaller subpopulation to the Endo2 poorly described signature that was detected primarily in lung and skin tumors [18] (Fig. 3F). Similar subpopulations of endothelial cells were also detected in SiNET2, highlighting the Endo2 signature as a recurrent pattern of subsets of SiNET endothelial cells (Fig. S3E-F).

Among B cells, two distinct populations were found in SiNET2 (Fig. 3G). Differential expression showed that the small B cell subpopulation is distinguished by upregulation of many cell cycle genes including canonical markers (MKI67, TOP2A, CDK1), reflecting proliferating B cells (Fig. 3H-I). Moreover, we noticed that, among all SiNET2 cells, the canonical cell cycle markers are expressed in B cells more than in all other cell types, including the malignant NE cells (Fig. 1C). This observation prompted us to turn our attention to cell cycle patterns across all cell types within the SiNETs.

Proliferation of NE and immune cells in SiNETs

Most SiNETs are low-grade tumors with a low mitotic index, but the exact identify of proliferating cells in SiNET is unknown. We used previously defined signatures of the cell cycle to identify all cycling cells. Notably, as we and others demonstrated previously [18, 22, 23], cell cycle involves the consistent upregulation of dozens of canonical genes and therefore the cycling cells can be robustly detected by scRNA-seq along with their phase along the cell cycle (Fig. S4).

Surprisingly, extremely few cycling cells were observed among the malignant NE cells (0.246% on average) (Fig. 4A). This fraction is considerably lower than what we detect in other cancers types [22], and as an extreme example from the same dataset, in the SiAdeno sample we find ∼13% of cycling epithelial cells using the same method (Fig. S4). Notably, the fraction of SiNET cycling cells was lower for NE cells than for all other cell types identified in the SiNET ecosystem (Fig. 4A; see also Fig. 1C,D). Relatively high fractions of cycling cells (>10%) were found only in epithelial or in B/plasma cells. Epithelial cells were only detected in two tumors, and of those only in one tumor they had a high fraction of cycling cells. In contrast, B/plasma cells were detected more commonly – in eight tumors – and in three of those they had a high fraction of cycling cells. The other cell types all had low-to-intermediate fractions of cycling cells across all SiNETs, but even those fractions were consistently higher than for the NE cells. Thus, at least in our SiNET cohort, proliferation is associated with the tumor microenvironment rather than with the NE cells. This result raises intriguing questions regarding the manner by which SiNETs grow, and the meaning of their mitotic index (see Discussion).

Figure 4.

Next, to ensure that this result is not unique to our cohort we reanalyzed bulk RNA-seq data of 81 SiNET samples [15]. We reasoned that if NE cells are a major source of proliferation, then cell cycle signatures should correlate with the expression of NE marker genes across bulk SiNET samples. Similarly, if other cell types, such as B cells or T cells, are more proliferative than the NE cells, then their markers should correlate with the cell cycle signature. As expected from the single cell analysis, we found that cell cycle signature significantly correlates positively with B cell markers, and to a more limited degree with T cell markers, but not with NE markers (Fig. 4B), supporting the broader relevance of our results.

B cell proliferation and MIF upregulation in the epithelial-like SiNET subtype

One potential explanation for the high proliferation of B cells is that our siNET samples may have included germinal centers (GC), in which B cells are expected to proliferate. To assess this possibility, we scored the B cells for a previously defined GC signature [24]. While two tumors with high B cell proliferation also had high GC scores, this was not the case for the third tumor (Fig. S5). Thus, inclusion of GCs may partially explain the unusual proliferation of B cells.

An alternative possibility is that B cells proliferate in response to particular signals in the SiNET microenvironment, for example, due to factors secreted by other cells. The proliferation of B/plasma cells was high only in three of the SiNETs and was absent or low in five other SiNETs (Fig. 4A). Intriguingly, all of the epithelial-like SiNETs had high B/plasma cell proliferation and all of the neuronal-like SiNETs (in which B/plasma cells were detected) had low proliferation of those cells. This perfect agreement with SiNET subtypes suggests that the unique features of the epithelial-like subtype may drive the proliferation of B/plasma cells. To identify such features, we defined the differential expression between the two subtypes for each cell type (Table S4).

This analysis highlighted Macrophage Migration Inhibitory Factor (MIF) as a prominent marker of the epithelial-like subtype that could potentially drive the proliferation of B/plasma cells. First, MIF was one of the top markers of the epithelial-like subtype in our analysis of NE cells (Fig. 2C), and a closer inspection reveals an extreme degree of differential expression, with almost no detected reads in the neuronal-like subtype and a consistently high expression in the epithelial-like subtype (Fig. 4C). Second, upregulation of MIF in the epithelial-like subtype was also detected all other cell types, although the effect was strongest in NE cells (Fig. 4C). Such consistent upregulation across six different cell types is not seen for any other genes, highlighting the unique upregulation of MIF in the ecosystem of the epithelial SiNETs, and supporting its potential causal effect in driving other phenotypes such as B/plasma cell proliferation. Third, MIF is an important cytokine that binds CD74 and was previously shown to regulate B cell proliferation by multiple studies [25, 26].

Putative progenitors in mixed tumors

Apart from our focus on SiNET, we were intrigued by the phenomena of MiNEN and aimed to profile such tumors by single cell RNA-seq. Due to their rarity and the difficulty of recognizing mixed tumors prior to their surgery, we were so far able to profile only one such tumor. This tumor was a Large Cell Neuroendocrine Carcinoma of the lung (LCNEC), mixed with a squamous cell carcinoma. While this tumor is highly different from SiNETs, we believe that this N-of-1 case study is of interest and raises important questions that might also be of relevance for SiNET and other NENs and hence we describe its analysis below.

We profiled the fresh LCNEC tumor sample by scRNA-seq using the 10x Chromium. Cells were first classified as malignant or non-malignant, based on inference of copy-number aberrations, as described previously [27] (Fig. 5A). Clustering of the non-malignant cells and annotation of the clusters based on standard markers separated them into fibroblasts, endothelial cells, macrophages, T cells and B cells (Fig. 5B).

Figure 5.

We found extensive diversity among the malignant cells, which highlighted several subpopulations of cells (Fig. 5B-C). These include the two malignant components expected based on the classification as a mixed tumor – neuroendocrine cells (e.g. expressing CHGB and SSTR2) and epithelial squamous cells (e.g. expressing KRT5 and KRT8). The epithelial cells further varied in their expression of a program that we previously termed “epithelial senescence” [28], which was indeed associated with lack of proliferation (Fig. 5D). Interestingly, we also found additional subpopulations of malignant cells with intermediate expression levels of both the neuroendocrine genes and the squamous genes (Fig. 5E). The “intermediate” cells included two clusters that we termed undifferentiated and progenitors, since they were distinguished by expression of developmental and stemness-related regulators (Fig. 5C). These include the homeobox transcription factors PBX1 and PROX1 [29, 30], The RNA-binding protein MSI2 [31], the WNT-related TCF4 [32, 33], and the NOTCH-related transcription factor HES6 [34] (Table S5). All four malignant populations were highly proliferative, except for the subset of squamous cells that express the epithelial senescence program (Fig. 5D). The progenitor cluster, with intermediate expression of the two lineage programs and upregulation of stemness factors, is consistent with the possibility that stem/progenitor cells may give rise to the two differentiated components that constitute the mixed phenotype. Additional MiNENs would need to be profiled at single cell resolution to explore the generality of this observation.

Discussion

While low-grade SiNETs are known to have a low mitotic index, any cancer is driven by cell proliferation such that malignant cells are expected to have higher cell cycle activity than non-malignant cells. Surprisingly, this does not seem to be the case in our analysis of SiNETs, where malignant cells appear to be less proliferative than multiple types of non-malignant cells. In such cases, it is conceivable that pathological evaluation of tumor proliferation (Ki67 counts) may be confounded by non-malignant cycling cells and therefore inflated. But a more fundamental question is how could those tumors initiate and progress when the NE cells appear to be largely non-proliferative?

We can envision multiple potential answers although further studies are needed to resolve these possibilities. First, since growth is determined by the balance between cellular proliferation, death and migration, SiNETs may still grow despite their very low proliferation through abnormally low death or abnormal patterns of migration. We are not aware of any evidence in support of this possibility but cannot exclude it. Second, proliferation of NE cells may be inhibited by prior treatments with Somatostatin analogues. Third, proliferation may be spatially or temporally restricted. Accordingly, higher proliferation may occur only in particular niches or at a particular time of tumor progression. This possibility would require an explanation for why our single cell analysis have consistently missed the relevant context of proliferation across the ten SiNETs.

A related possibility is that proliferation may be restricted to a hidden stem/progenitor population that we have not captured, either because it is rare or because it may be sensitive to our experimental workflows. Such proliferative stem/progenitor cells were described in other cancer types and are consistent with the cancer stem cell hypothesis. For example, we previously described the cellular hierarchy of oligodendroglioma, in which proliferation is restricted to cells resembling neural progenitor, while the tumor is dominated by more differentiated cells resembling oligodendrocytes and astrocytes [27]. A second, and more relevant example, is presented here with the N-of-1 analysis of a mixed LCNEC. In this case, high proliferation is seen across all malignant components, but the mere presence of a proliferative stem/progenitor population suggests that bipotent progenitors may explain the diversity within mixed LCNECs and possibly even the wider phenomena of MiNENs. It is tempting to speculate that such progenitors may also exist in other neuroendocrine tumors (e.g. SiNETs).

SiNETs deviate from classical oncogenesis not only by their low malignant cell proliferation, but also by their paucity of mutations and apparent driver events. Accordingly, previous work proposed that SiNET may be driven by epigenetic aberrations in NE cells or by effects of the TME. Our observation of minimal proliferation in NE cells supports the latter possibility of TME oncogenic effects. In particular, SiNETs of the epithelial-like subtype are associated with extremely high levels of MIF, an important and pleiotropic cytokine that may influence the TME in multiple important ways and possibly create a unique inflammatory TME. One potential effect that we observe is high proliferation of B or plasma cells. Notably, the three SiNETs of the epithelial-like subtype are not associated with higher numbers of B/plasma cells compared to other SiNETs, but rather only with increased cell cycle of B/plasma cells, suggesting that their proliferation is balanced by high death. The implications of high B/plasma cell turnover, and of other downstream effects of high MIF expression, are unclear, but raise the possibility that MIF-CD74 interaction may constitute a relevant target for the epithelial-like SiNET subtype. Notably, MIF-CD74 inhibitors have already been developed and initially tested for the treatment of solid tumors [35-37].

In summary, single cell profiling of SiNETs revealed two distinct subtypes that differ in NE profiles and in the proliferation of B-cells, possibly justifying different therapeutic strategies (e.g. MIF inhibition). The two subtypes share a minimal proliferation of NE cells, which raises questions about the initiation and growth of SiNETs and the potential role of TME cells, in line with the limited mutations observed in SiNETs. Our analysis of a mixed tumor further suggests the potential presence of bipotent progenitors that may account for mixed phenotypes. Taken together, these results improve our understanding of neuroendocrine neoplasms while highlighting their complex biology.

Acknowledgements

This work was supported by grants from the Neuroendocrine Tumor Research Foundation (to I.T. and A.T.). I.T. is further supported by the Zuckerman STEM Leadership Program, the Mexican Friends New Generation, the Benoziyo Endowment Fund and E. Harari. I.T. is the incumbent of the Dr. Celia Zwillenberg-Fridman and Dr. Lutz Zwillenberg Career Development Chair.

Data availability

Data will be available through GEO.

Competing interests

I.T. is advisory board member of Immunitas Therapeutics.

Methods

Human samples

The study was approved by the Instituntional Helsinkee Committee at Sheba Medical Center (SMC-18-5674). For tissue preparation for Immunohistochemistry, flash frozen SiNETs were stored at −80°C until cutting. Blocks were prepared in OCT followed by freezing and sectioning. Frozen tissue was sectioned using a −20°C temperature on a cryostat (Leica CM3050 S) onto microscopic slides (Thermo, Superfrost Plus).

Tumor dissociation and library preparation

Fresh tumors were collected directly from the operating room at the time of surgery. Tumors were washed in ice-cold HBSS, minced using a pair of Noyes spring scissors (Fine Science Tools) and then enzymatically dissociated using a tumor dissociation kit (Miltenyi Biotec) in a GentleMACS Octo-dissociator on low speed settings (Filbin et al., 2018; Patel et al., 2014; Tirosh et al., 2016b). Single cell suspension was obtained by passing dissociated slurry through a 60-100 micron cell strainer (Miltenyi) and the filtrate was subjected to RBC removal using a RBC lysis buffer (Roche). Dead cell removal (Miltenyi) was performed on cell suspension following manufacturer’s protocol and cells were assessed for viability and density. 5 µl of Trypan blue (Thermo Fisher Scientific) was mixed with 5 µl of the sample and loaded onto chip disposable automated hemocytometer (Countess II). Cell concentration was adjusted such that a total of 8000-10,000 cells were loaded onto each channel of the 10x Genomics Single-Cell Chromium Controller.

For frozen tumors, a similar dissociation procedure was performed with the following modifications. Nuclei isolation was achieved using either ST-based buffers as done before (Slyper et al) or by EZ-lysis method (Sigma). Briefly, tissue samples were cut into pieces <0.5 cm using spring noyes scissors and then homogenized in ice-cold buffer (ST based or EZ lysis) using a glass Dounce tissue grinder (Sigma). Nuclei were centrifuged at 500g, with low acceleration/deceleration for 5 min at 4 °C, washed with 5 ml ice-cold ST/EZ lysis buffer and incubated on ice for 5 min. A suspension buffer containing PBS, BSA and RNAse inhibitor was used to wash and resuspend the nuclei pellet. Nuclei suspension was filtered through a 40 μm cell strainer (Miltenyi) and counted. A final concentration of 1,000 nuclei per µl was used for loading on a 10x channel with a target of 10,000 nuclei per channel in the 10x controller.

The Chromium Next GEM Single Cell 3′ GEM, Library & Gel Bead Kit v3.1, Chromium Single Cell 3′ Feature Barcode Library Kit, Chromium Next GEM Chip G, and 10x Chromium Controller (10x Genomics) was used for generation of libraries of fresh tumors/single cells and the Chromium Next GEM Single 5′ GEM was used for frozen tumors/single-nuclei.

As per the standardized protocol of creating 10x libraries, single cells, reagents and single gel beads containing barcoded oligonucleotides were emulsified and encapsulated into nanoliter-sized droplets followed by reverse transcription (RT). Following manufacturer’s recommendations RT samples were subjected to cDNA amplification, fragmentation, adapter and sample index attachment. Libraries from two 10x channels were pooled together and sequenced on one lane of an Illumina NovaSeq-6000, using an sp-100 kit or 4 were pooled together on two lanes with an S1-100 kit, with paired end reads as follows: read 1, 26 nt; read 2, 55 nt; index 1, 8 nt; index 2, 0 nt.

Sample Normalization, Filtering and annotations

Four samples (siAdeno, SiNETs 1-3) underwent single nuclei sequencing, while the remaining samples (SiNETs 4-10) were subjected to scRNA-seq. To facilitate comparative analysis, we converted counts of unique molecular identifier (UMI) counts to transcripts per million (TPM) values, which were divided by 10 since the actual complexity of cells is assumed to be in the realm of ∼100,000 transcripts and not 1 million as implied by the TPM measures. The resulting values were added 1 and log2-transformed. Finally, to derive relative expression we centered the value of each gene in each tumor, except in analysis where we were interested in absolute expression levels.

To ensure data quality, we excluded low-quality cells based on the number of detected genes. For the nuclei-seq platform, we used a threshold of 1000 detected genes, while for the scRNA-seq platform, we used a threshold of 700 in order to retain lymphocytes that tend to have low numbers of detected genes. Additionally, we performed most analysis with a filtered set of genes, retaining only the genes in which the logged row-means per gene across all cells passed 4.5. After filtering, the number of cells analyzed per sample ranged from 287 to 10,802.

Neuroendocrine cell signatures and their classification

To identify NE-specific genes, we performed a differential expression analysis between NE cells and a reference group consisting of Macrophages, Fibroblasts, and Endothelial cells. For each sample, we sampled 50 NE cells along with 50 cells from the reference group. Differentially expressed (DE) genes were defined as those with fold change greater than 8 and a p-value lower than 0.05. In samples that had insufficient reference cells (SiNET7 and SiNET8), we sampled the reference cells from a different sample that exhibited the highest correlation with the given sample. Samples that did not contain NE cells (SiNET3 and siAdeno) were excluded from this analysis. We then classified the DE genes into three groups: Common genes, Subset-specific genes, and Sample-specific genes. Common genes were defined as those differentially expressed in at least 6 samples, Subset-specific genes were those differentially expressed in 2-5 samples, and Sample-specific genes were differentially expressed only in one sample.

Clustering of Subset-specific NE genes into nine patterns

We divided the subset-specific NE genes based on their expression pattern across the 8 tumors that contained enough NE cells. To this end, we first considered all potential binary patterns, in which a gene is considered as either expressed (1) or not expressed (0) in each tumor. This defined 2⁸ = 256 theoretical patterns. Next, for each gene, we calculated the correlation between its relative expression across the eight samples, and each of the 256 binary patterns; the gene was then assigned to the pattern with highest correlation. This uncovered 9 patterns that were each assigned to >5 genes (Fig.2C, Table S3).

Cell Cycle Analysis

To investigate the cell cycle dynamics within the dataset, we utilized a canonical cell cycle gene-set [27]. Each cell in every sample was scored based on the expression of these cell cycle genes. To identify cycling cells, we established a threshold based on the distribution of scores within each sample. We chose a relatively lenient threshold (log2 fold change of 1.5), and further verified that our main claims are qualitatively maintained when using a more strict threshold of 2-fold. In particular, NE cells had an extremely low fraction of cycling cells (0.127%), while B and plasma cells had a high fraction of cycling cells in epithelial-like SiNETs (32.25%) but not in Neuronal-like SiNETs (0.178%).

Within Cell-Type analysis

In each tumor, we examined clusters that were annotated as the same cell type but were not grouped together in the initial clustering analysis. For each such pair of clusters, we defined differential expression (with log-ratio of 1.5 and a p-value of 0.05), and analyzed the expression of previously defined meta-programs associated with the corresponding cell type.

Comparing Gene Expression between SiNET Subtypes across Cell Types

We performed a comparative gene expression analysis, examining epithelial-like and neuronal-like subtypes of SiNETs. The significance threshold for gene selection was established at an uncorrected p-value < 0.05, along with a log2 fold change > 2. Given the limited sample size. Additionally, we extended this analysis to bulk data, where each sample was assigned to the subtype for which the signature score was higher.

Cell cycle correlations with cell type signatures in bulk SiNET samples

From GSE98894, we used the 81 SiNET samples. Counts were converted to TPM, and log2-transformed with an offset of 1. We then averaged the marker genes of Epithelial, T cells, B cells, NE cells and cell cycle to create scores of their expression. Few outlier samples with very low NE score were excluded from further analysis. Finally, we calculated Pearson correlation between the cell-type scores and cell cycle score.

Supplementary figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.