Abstract

   Neuroendocrine tumors (NETs) are rare cancers that most often arise in
   the gastrointestinal tract and pancreas. The fundamental mechanisms
   driving gastroenteropancreatic (GEP)–NET growth remain incompletely
   elucidated; however, the heterogeneous clinical behavior of GEP-NETs
   suggests that both cellular lineage dynamics and tumor microenvironment
   influence tumor pathophysiology. Here, we investigated the single-cell
   transcriptomes of tumor and immune cells from patients with
   gastroenteropancreatic NETs. Malignant GEP-NET cells expressed genes
   and regulons associated with normal, gastrointestinal endocrine cell
   differentiation, and fate determination stages. Tumor and lymphoid
   compartments sparsely expressed immunosuppressive targets commonly
   investigated in clinical trials, such as the programmed cell death
   protein–1/programmed death ligand–1 axis. However, infiltrating myeloid
   cell types within both primary and metastatic GEP-NETs were enriched
   for genes encoding other immune checkpoints, including VSIR (VISTA),
   HAVCR2 (TIM3), LGALS9 (Gal-9), and SIGLEC10. Our findings highlight the
   transcriptomic heterogeneity that distinguishes the cellular landscapes
   of GEP-NET anatomic subtypes and reveal potential avenues for future
   precision medicine therapeutics.
     __________________________________________________________________

   Neuroendocrine tumors bear hallmarks of early gastrointestinal
   development and an immunosuppressive myeloid microenvironment.

INTRODUCTION

   Neuroendocrine tumors (NETs) are rare cancers of the diffuse
   neuroendocrine system, and while they represent less than 1% of all
   newly diagnosed malignancies in the United States per year, the
   population prevalence is estimated to exceed 100,000 individuals
   ([92]1–[93]3). The diagnosed incidence of GEP-NETs has risen markedly
   over the past four decades, a persistent trend ascribed partially to
   improved tumor classification and diagnostic technologies ([94]1,
   [95]2, [96]4). Approximately two-thirds of all NETs are
   gastroenteropancreatic (GEP) in origin, predominantly in the pancreas
   (12 to 30% of GEP-NETs), stomach, and small intestine (31 to 50% of
   GEP-NETs) of older adults ([97]3, [98]5). While patients with localized
   disease can be cured with surgical resection, these tumors are often
   identified at a late clinical stage. Despite rapid advancements in NET
   research and clinical investigation, treatment options for such
   patients remains limited ([99]5).

   A limiting factor in developing novel therapies for NETs is that the
   fundamental biology of GEP-NETs remains incompletely understood.
   Previous investigations have been challenged by limited tissue
   availability and a scarcity of representative in vitro and in vivo
   experimental models. Characteristic genomic mutations in these tumors
   most often occur in tumor suppressor genes associated with hereditary
   endocrine cancer syndromes, including multiple endocrine neoplasia type
   1 (MEN1), von Hippel-Lindau disease (VHL), or tuberous sclerosis
   (TSC1/2) ([100]6). Chromosomal loss, telomere alterations, and
   epigenetic modifications have emerged as potential regulators of
   GEP-NET development but are not sufficiently targeted by currently
   available treatments ([101]7, [102]8). The inter- and intratumoral
   heterogeneity of GEP-NETs at different anatomical sites (e.g.,
   pancreas, small intestine, and stomach) poses an additional challenge
   to understanding the oncogenic mechanisms driving tumor formation in
   these disparate environments.

   High-resolution delineation of the tumor and immune niches of GEP-NETs
   has the potential to elucidate previously unknown, actionable targets
   for future precision medicine endeavors, expanding the limited menu of
   therapeutic options available to patients with nonresectable tumors.
   Previous investigations to this end have used immunohistochemical or
   bulk RNA sequencing approaches, but these methodologies are
   insufficient to identify rare or unknown cell states that contribute to
   tumor biology and progression. Single-cell RNA sequencing (scRNA-seq)
   of other gastrointestinal solid tumors has enabled high-resolution
   characterization of the cell subpopulations within the cancer
   microenvironment, thereby overcoming the challenges presented by older
   technologies ([103]9–[104]12). We therefore hypothesized that
   dissecting the transcriptomic landscape of GEP-NETs at single-cell
   resolution would elucidate shared and distinct cellular phenotypes and
   genetic programming among these heterogeneous tumors.

RESULTS

Characterizing the cellular landscape of GEP-NETs

   We collected fresh surgical resection specimens from eight patients
   with well-differentiated NETs, comprising four pancreatic NETs (pNETs),
   three small intestinal NETs (siNETs), and one gastric NET (gNET). The
   siNET and gNET samples included tissue from regional extension of the
   primary tumor. Resection tissue was dissociated and divided in two for
   parallel processing; one-half of each cell suspension sorted for
   viability using flow cytometry for scRNA-seq using the 10x Genomics
   platform, while the other was flow sorted for both CD45 positivity and
   viability before undergoing 10x processing ([105]Fig. 1A and Methods).
   Of the eight samples in our cohort, five had received no previous
   therapy, while the remaining three tumors had been treated with the
   somatostatin analog lanreotide. The gNET had been treated with both
   lanreotide and everolimus and chemotherapy with capecitabine and
   temozolomide ([106]Fig. 1B). Our cohort also included one
   primary-metastasis pair from the same patient (sinet2 and sinet3;
   [107]Fig. 1B).

Fig. 1. Characterizing the cellular landscape of pancreatic, small
intestinal, and gNETs.

   [108]Fig. 1.
   [109]Open in a new tab

   (A) NET scRNA-seq workflow. Resected tumor fragments were dissociated
   into single-cell suspensions, which were then divided for parallel
   sequencing of (i) nonsorted single cells and (ii)
   fluorescence-activated cell sorting (FACS) sorted, CD45^+ immune cells.
   All samples were processed using the 10x Genomics platform. (B) Patient
   demographics and treatment summary at time of tumor resection. Site:
   SI, small intestine; SI/liver: small intestinal primary tumor with
   liver metastasis. A bar underneath samples indicates that they
   originate from the same patient. (C) UMAP of malignant and nonmalignant
   cells from all tumors, colored by general cell type. (D) UMAP of
   malignant and nonmalignant cells from all tumors, colored by sample of
   origin (left). (E) Cellular composition of each tumor, colored by
   general cell type.

   After preprocessing and quality control, our single-cell atlas
   contained high-quality transcriptomes from 24,048 cells spanning both
   tumor and immune compartments ([110]Fig. 1C and Methods).
   Across-dataset integration of all cells (CD45^+ sorted and non–CD45^+
   sorted) was performed using the Harmony algorithm to mitigate batch
   effects stemming from technical or biological variances between
   individual patients, tissue collection, and sample processing
   ([111]Fig. 1, C and D, and Methods). Using differential gene expression
   analysis, we identified cluster-specific marker genes that enabled us
   to subcategorize the lymphoid and myeloid immune compartments
   ([112]Fig. 1C and fig. S1A). Tumor cells were identified as cells
   differentially expressing NET-related marker genes and demonstrating
   inferred copy number variation compared to reference normal cells (fig.
   S1C and Methods).

   We then compared the microenvironmental composition of GEP-NETs from
   distinct sites of origin ([113]Fig. 1E). In general, myeloid cells
   including macrophages, monocytes, and dendritic cells (DCs) represented
   the largest immune subpopulation within pNETs, while siNETs and the
   gNET contained a relatively higher proportion of T cells and natural
   killer (NK) cells ([114]Fig. 1E). Cell type proportions varied between
   GEP-NETs from the same site of origin, reflecting the intrinsic
   heterogeneity of these tumors.

Differential expression analysis and gene regulatory network inference reveal
multilineage profiles of pNETs

   We next interrogated the transcriptional landscapes of GEP-NET cells
   from different anatomical sites of origin ([115]Fig. 2A). All tumors
   expressed known NET markers such as chromogranin genes CHGA and CHGB
   ([116]Fig. 2B). While most tumors expressed somatostatin receptor 2
   (SSTR2), no GEP-NET expressed SSTR5 ([117]Fig. 2B). We also
   interrogated commonly mutated genes in pNETs to assess their expression
   levels across GEP-NETs. ATRX, a somatic mutation identified in pNETs,
   was heterogeneously expressed across GEP-NET subtypes; DAXX,
   conversely, was either minimally or not expressed by any tumor
   ([118]Fig. 2B). Mutations of ATRX and DAXX in pNETs are often
   inactivating or missense in nature, which could explain the paucity of
   expression we identified here ([119]13).

Fig. 2. GEP-NET cells do not express common immune checkpoint targets.

   [120]Fig. 2.
   [121]Open in a new tab

   (A) UMAP plot of all malignant cells from all GEP-NET samples, colored
   by sample of origin. (B) Violin plots of genes commonly expressed in
   GEP-NETs, grouped and colored by sample. (C) Heatmap of scaled,
   normalized average expression of top 50 differentially expressed genes
   between pancreatic (pn) and small intestinal (si) NETs. (D) Heatmap of
   scaled expression of the top five differentially expressed regulons per
   individual tumor. Blue-labeled genes indicate involvement in
   enteroendocrine differentiation or regulation of intestinal stem cell
   fate. Red-labeled genes indicate TFs with known involvement in
   endocrine pancreas development and pNET biology. (E) UMAP plots of
   malignant cells from all pNETs, colored by tumor of origin (top),
   followed by heatmaps in UMAP space of VISION signature scores for cell
   type–specific genes for the major cells in the normal human endocrine
   pancreas (bottom). (F) Dotplot of percent expression of common,
   targetable immune checkpoint and evasion genes, grouped by sample. (G)
   Representative immunohistochemical staining of CD8, PD-1, and PD-L1 two
   siNET samples. PD-1 and PD-L1 expression was low to nonexistent across
   all stained samples, consistent with low transcriptional expression by
   scRNA-seq. Scale bars, 100 μm. DAPI, 4′,6-diamidino-2-phenylindole.

   To compare the transcriptional profiles of pNETs and siNETs, we then
   performed differential gene expression on all tumor cells from those
   anatomical sites. Primary and metastatic siNETs exhibited relatively
   homogeneous transcriptional landscapes compared to pNETs and featured
   genes involved in monoamine and hormone synthesis (e.g., TAC1 and TPH1)
   and neuroendocrine transport and vesicle release (e.g., SYT13 and
   SLC18A1) ([122]Fig. 2C), concordant with established profiles of
   carcinoid tumors. The pNETs in our cohort displayed marked intertumoral
   heterogeneity, with few differentially expressed genes common to all
   four tumors ([123]Fig. 2C). However, all pNETs expressed genes involved
   in pancreas exocrine (e.g., CCK and SPINK1) and endocrine function
   (e.g., INS, VIP, GCG, and PPY) ([124]Fig. 2C). The neuroendocrine genes
   expressed by each tumor corresponded closely to its clinical profile;
   for example, pnet2 demonstrated the highest expression of the
   insulin-encoding gene INS, matching its pathological characterization
   as a well-differentiated insulinoma ([125]Fig. 2C).

   Next, we sought to identify regulatory elements underlying the
   heterogeneous molecular profiles of GEP-NETs from different anatomical
   sites. To this end, we performed single-cell regulatory network
   inference and clustering (SCENIC) analysis on malignant cells from each
   tumor. We identified the top regulons for each tumor by the area under
   the curve (AUC) value as calculated by the standard pySCENIC workflow
   ([126]Fig. 2D and Methods). Hierarchical clustering of NETs by regulon
   activity produced two major subgroups that did not segregate solely by
   anatomical location; the first contained most pNETs and the gNET
   sample, while the other comprised all siNETs and pnet1. Within our pNET
   samples, one tumor (pnet2) displayed differential activity of the
   transcription factors (TFs) PDX1, PAX6, MAFA, NKX6-1, and RXRG
   ([127]Fig. 2D). These TFs are known regulators of beta cell fate in
   healthy human pancreatic islets, consistent with the clinical diagnosis
   of insulinoma ([128]14–[129]17). Another pNET (pnet1) displayed
   differential activity of the ARX TF, which primarily directs pancreas
   alpha cell development but also has known roles in establishing
   pancreatic gamma cell fate in healthy islets of Langerhans
   ([130]Fig. 2D) ([131]18, [132]19). These data are consistent with a
   previous investigation of enhancer profiles of nonfunctional pNETs,
   which found that those tumors could be stratified into two major groups
   by PDX1 and ARX activity ([133]19). However, two of our pNETs did not
   express either PDX1 or ARX. Pnet4 displayed differential activity of
   TCF12, which has been implicated in pancreatic and gastrointestinal
   cancer proliferation and invasion, while pnet3’s top regulon, RARB, was
   previously found to direct differentiation of embryonic stem cells to
   pancreatic endocrine cells in vitro ([134]20, [135]21). Few common TFs
   were identified across NET samples, revealing the heterogeneity within
   and between GEP-NETs from different anatomical sites.

   Within our siNET samples, we identified differential activity of TFs
   associated with regulation of enteroendocrine cell lineages (NEUROD1
   and FOXA1) and with the maintenance of intestinal stem cell fate
   (KLF5), concordant with a leading hypothesis in the field of NET
   biology that siNETs develop from enterochromaffin cells in the gut
   ([136]Fig. 2D) ([137]22–[138]24). The lone siNET metastasis sample
   (sinet3) displayed highly differential activity of homeobox (HOX) TFs,
   including HOXB6 and HOXB13 ([139]Fig. 2D). Previous analysis of primary
   human fetal small intestinal tissue revealed that HOXB gene expression
   is strongest early in duodenal development but quickly diminishes over
   time ([140]25). The HOXB family was not represented amongst the top
   regulons for either primary siNET in our cohort, suggesting that the
   siNET metastasis may have undergone dedifferentiation to resemble an
   earlier stage of small intestinal fate determination ([141]Fig. 2D).

   To determine the transcriptional similarity of our four pNET tumors to
   healthy, mature pancreatic islet cells, we scored malignant cells from
   each tumor with curated gene signatures specific to alpha, beta, gamma,
   and delta cells in mature human pancreatic islets ([142]Fig. 2E)
   ([143]26). Tumor cells from pnet1 displayed robust expression of alpha
   cell and gamma cell gene signatures, with lower expression of the delta
   cell signature ([144]Fig. 2E); these findings are consistent with the
   robust expression of gamma and delta cell hormones PPY and VIP and
   elevated ARX TF activity identified in our previous analyses.
   Similarly, as predicted by the results of our differential gene
   expression and gene regulatory network (GRN) analysis, pnet2 was the
   tumor with the highest expression score for the beta cell signature
   ([145]Fig. 2E). Pnet3 and pnet4 scored most highly for the alpha cell
   signature with lower expression of delta and gamma cell scores
   ([146]Fig. 2E). Together with the results of our regulatory network
   analysis, these data demonstrate tumor cell resemblance to earlier
   stages of gastrointestinal development and the potential for lineage
   plasticity in these rare tumors.

GEP-NET cells do not express classical immune checkpoint blockade targets

   The success and ongoing studies of immunotherapy in other solid tumor
   types have motivated a strong interest in applying this therapeutic
   approach to the treatment of GEP-NETs. We therefore next sought to
   quantify the expression of classical, tumor cell–anchored immunotherapy
   targets across the GEP-NETs ([147]Fig. 2F). Programmed death ligand–2
   (PDL2), which binds programmed cell death protein 1 (PD-1) on the
   surface of T cells to suppress antitumoral immunity, was not expressed
   in any GEP-NETs in our cohort. PDL1, which also binds to PD-1 to
   negatively regulate T cell immunity, was expressed at low levels in
   fewer than half of tumors. Although roughly 40% of gNET cells expressed
   HLA-G, a recently described immunomodulatory ligand, all other tumors
   rarely expressed this gene; anatomical location, pathological grade,
   and treatment history may therefore drive this differential expression
   ([148]27). LGALS9 (also known as Gal-9), VSIG, and VSIR, other
   regulators of T cell anergy and exhaustion, were also sparsely
   expressed but were identified in a greater proportion of tumors than
   PDL1.

   Next, we sought to validate the low expression of the PD-1/PD-L1
   immunosuppressive axis in GEP-NETs at the protein level. To do this, we
   performed multiplex immunofluorescence staining of PD-1 and PD-L1 in
   formalin-fixed, paraffin-embedded (FFPE) tissue from small intestinal
   and pNETs in our cohort using the ImmunoProfile platform ([149]28).
   Across all stained samples, we identified very few PD-1– or
   PD-L1–positive cells in the tumor microenvironment and found no
   positive cells in the tumor compartment ([150]Fig. 2G and fig. S2, A to
   C). Together, our observations can inform future clinical investigation
   into the use of anti–PD-1/PD-L1 immunotherapies versus emerging
   immunotherapy options in patients with well-differentiated GEP-NETs.

   Given the demonstrated efficacy of the mammalian target of rapamycin
   (mTOR) inhibitor everolimus in the treatment of patients with advanced
   GEP-NETs, we also investigated the expression of mTOR pathway genes in
   our cohort ([151]29). We scored all malignant cells with the
   “PI3K/AKT/mTOR Signaling” gene signature from the HALLMARK collection
   in Molecular Signatures Database (MSigDb; fig. S2, D and E) and
   performed pairwise score comparisons to establish the significance of
   expression differences (Wilcoxon two-sided test, P < 0.05). The siNETs
   and gNET in our cohort scored significantly higher for the
   PI3K/AKT/mTOR Signaling signature compared to the pNETs (fig. S2D). The
   metastatic tumor sinet3 had a lower signature score than its primary
   counterpart (sinet2). Functional follow-up will be required to
   determine the correlation between gene signature score and everolimus
   response across GEP-NET subtypes, but our findings demonstrate
   heterogeneity in mTOR signaling across GEP-NETs from different sites of
   origin.

Profiling lymphoid cell heterogeneity in GEP-NETs

   To better understand the immune microenvironment of GEP-NETs, we next
   analyzed the lymphoid compartment of our dataset. We selected all T and
   NK cells from the full scRNA-seq atlas using canonical marker genes
   (fig. S1A and Methods) for de novo clustering and differential gene
   expression analysis (fig. S3, A and B). We identified diverse T and NK
   cell populations including NK cells (FGFBP2^+ and FGFBP2^−), resting
   and effector CD4^+ T cells, regulatory T cells (T[regs]), CD8^+ T
   cells, T helper 17 cells (T[H]17), and mucosal-associated invariant T
   cells (MAITs) ([152]Fig. 3, A and B). These T and NK cell types were
   identified in differing proportions across GEP-NET site of origin, with
   the gNET containing the highest proportion of T[regs] across samples
   ([153]Fig. 3A). We also classified B cells from our original dataset,
   the majority of which originated in either the gNET or siNETs (fig.
   S3C). Of the six B cell clusters identified in our analysis, two
   contained IL4R-expressing naïve B cells, while the other four comprised
   CD27-expressing memory B cells (fig. S3, D and E).

Fig. 3. Classifying T cell and NK cell populations within GEP-NETs.

   [154]Fig. 3.
   [155]Open in a new tab

   (A) UMAP (top) of T cells and NK cells from all tumors, colored and
   labeled by lymphoid subtype. Barplot (bottom) displays T/NK cell type
   proportion by sample. (B) Heatmap of scaled, normalized expression of
   top 15 differentially expressed genes per T/NK cell subtype. Specific
   subtype-determining genes are annotated at right. (C) UMAP of CD8^+ T
   cells from all tumors, colored by CD8^+ T cell subtype (left) and NET
   origin site (right). (D) Expression heatmaps in UMAP space of naïve
   (TCF7 and CCR7), memory (IL7R), and exhausted (TOX) CD8 T cell markers
   as well as coinhibitory receptor TIM3 in CD8^+ T cells. (E) Dotplots of
   percent of T cell subtypes expressing known immunosuppressive genes,
   grouped and colored by NET site.

   The higher expression of immune checkpoint genes in 4-1BB-hi CD8^+ T
   cells compared to other CD8^+ T cells across all GEP-NETs motivated us
   to investigate the exhaustion phenotypes present in infiltrating CD8^+
   T cells. For this analysis, we included all cell types with nonzero
   expression of CD8A or CD8B: 4-1BB-hi CD8^+ T cells, 4-1BB-lo CD8^+ T
   cells, and MAITs. We first isolated all CD8^+ cells from our T and NK
   cell dataset, which yielded two major subclusters ([156]Fig. 3C). One
   cluster contained mostly 4-1BB-hi CD8^+ T cells from the gNET sample,
   while the other cluster represented a mixture of all three CD8^+ T cell
   types from all GEP-NET sites of origin. To further delineate CD8^+ T
   cell subtypes, we interrogated the transcription expression patterns of
   genes previously identified in stem-like T cells (TCF7), memory T cells
   (IL7R and CCR7), activated T cells (CD27), and exhausted T cells (TOX
   and HAVCR2/TIM3) ([157]Fig. 3D). Expression levels of TCF7 and IL7R
   were highest in MAIT and 4-1BB-lo CD8^+ T cells, while TOX and TIM3
   displayed the exact opposite expression pattern. These findings are
   consistent with previous single-cell analyses in renal cell carcinoma,
   which identified high expression of exhaustion markers among 4-1BB-hi
   CD8^+ T cells ([158]30).

   To evaluate the potential therapeutic targets in the GEP-NET milieu
   using current immunotherapies, we quantified the expression of known,
   lymphoid-intrinsic immune checkpoint ligands on the T and NK cells in
   our atlas. CTLA4 and TIGIT were expressed robustly and frequently in
   T[regs] across all GEP-NET subtypes (i.e., pNETs, siNETs, and gNETs)
   ([159]Fig. 3E). LAG3, HAVCR2 (encoding TIM3), and PDCD1 (encoding PD-1)
   were sparsely expressed in T and NK cells from pNETs and siNETs. The
   lymphoid compartment of our one gNET sample demonstrated the highest
   expression of all queried immunosuppressive markers, with up-regulation
   in the 4-1BB-hi CD8^+ T cell, T[regs], FGFBP2^− NK cell, and cycling
   populations ([160]Fig. 3E). The gNET sample had the highest World
   Health Organization (WHO) Grade [grade 3 (G3)] of all analyzed samples
   and was the only tumor to receive treatment with everolimus and
   chemotherapy in addition to a somatostatin analog ([161]Fig. 1B), which
   may influence the up-regulation of immune checkpoint expression in its
   T and NK cells. Determining the roles of grade and treatment history in
   this phenomenon will necessitate expanded serial biopsy cohorts for
   future analysis.

Myeloid cells in GEP-NETs express targetable immunosuppressive ligands

   Given the limited expression of classical immune checkpoint ligands and
   receptors on lymphoid cells in the GEP-NET microenvironment, we next
   sought to characterize the myeloid cells within our cohort. We
   subset all myeloid cells (n = 4820 cells) from our original dataset
   using known marker genes and reclustered them to identify monocytes
   (CD16^+ and CD16^−), DCs (CLEC10A-high DCs, CLEC9A-high DCs, and
   CD207^+ DCs), and tumor-associated macrophages (FOLR2-high TAMs and
   SELENOP-high TAMs) ([162]Fig. 4A and fig. S4, A and B). The proportion
   of myeloid cell types was similar across all tumors, and we did not
   identify any patient-specific or GEP-NET subtype–specific myeloid
   populations ([163]Fig. 4A).

Fig. 4. GEP-NET–intrinsic myeloid cells express targetable immunosuppressive
ligands.

   [164]Fig. 4.
   [165]Open in a new tab

   (A) UMAP of myeloid cells (top) from all tumors, colored and labeled by
   myeloid cell subtype. Barplot (bottom) displays myeloid cell type
   proportion by sample. (B) Dotplot of percentage of each myeloid subtype
   expressing immunosuppressive ligand genes. (C) Two-dimensional
   butterfly plot visualizations of relative metamodule scores of top
   MSigDB Hallmark Pathways (TNF-α, “TNFA Signaling Via NFKB;” IFN-γ,
   “Interferon Gamma Response;” IL-2/STAT5, “IL2 STAT5 Signaling;” and
   complement = “Complement”) in all GEP-NET–associated macrophages (left)
   and separated by TAM subtype (right). (D) Violin plot of VISION scores
   of M1 macrophage– and M2 macrophage–associated gene signatures, grouped
   by TAM subtype. (E) Heatmaps of scaled, normalized, and averaged
   expression of top differentially expressed genes between GEP-NET sites
   within each TAM subtype (left, FOLR2-hi TAMs and right, SELENOP-hi
   TAMs). (F) Graphical summary of differential chemokine expression
   patterns between TAMs from GEP-NETs at different anatomical sites.

   To determine whether myeloid cells in GEP-NETs expressed targetable
   immune checkpoint ligands, we then quantified the expression of
   immunosuppressive genes across our nine myeloid clusters. As in the
   tumor compartment, PDL1 and PDL2 were expressed in fewer than 20% of
   any given myeloid subtype ([166]Fig. 4B). VSIG4, an inhibitor of both
   cytotoxic T cell and proinflammatory macrophage activation, was
   expressed most strongly in FOLR2-hi TAMs ([167]31, [168]32). Both TIM3
   and its binding partner Galectin-9 were more broadly expressed across
   both TAM and DC subtypes, suggesting an ability of these cells to
   suppress the activity of cytotoxic T cells. Last, the immunosuppressive
   ligands VSIR and SIGLEC10, which regulate antitumoral macrophage
   activity, displayed robust expression across the myeloid compartment;
   VSIR specifically was expressed in over 60% CD16^+ monocytes and in
   more than half of SELENOP-hi and FOLR2-hi TAMs ([169]33).

   To better define the diversity of tumor-infiltrating macrophage
   phenotypes in GEP-NETs, we isolated the four TAM clusters and performed
   single-cell pathway enrichment analysis using the HALLMARK genesets
   from the Molecular Signature Database (Methods). The most significantly
   enriched pathways included (i) tumor necrosis factor–α (TNF-α)
   signaling via nuclear factor κB, (ii) interferon-γ (IFN-γ) response,
   (iii) interleukin-2 (IL-2)/signal transducers and activators of
   transcription 5 (STAT5) signaling, and (iv) genes involved tumor growth
   factor–β signaling. Using the signature scoring protocol outlined in
   Neftel et al. ([170]34), we then scored each TAM for expression of each
   of the four geneset signatures, also referred to as metamodules. We
   visualized the results of this scoring analysis as a two-dimensional
   butterfly plot wherein each quadrant represents a specific metamodule
   ([171]Fig. 4C, left). FOLR2-hi TAMs scored highly across all four
   metamodules, while SELENOP-hi TAMs predominantly expressed TNF-α
   signaling and IFN-γ response genes ([172]Fig. 4C, center). Grouping
   TAMs by tumor site of origin revealed that while TAMs from pNETs span
   all four metamodules, TAMs from siNETs and the gNET predominantly
   differentially score for the IFN-γ response signature ([173]Fig. 4C,
   right).

   Next, we assessed whether GEP-NET TAMs could be defined according to
   the traditional, albeit debated, M1 versus M2 macrophage binary
   classification system ([174]35, [175]36). M1 macrophages are generally
   considered proinflammatory and antitumoral, whereas M2 macrophages
   display immunosuppressive activity. To this end, we scored all TAMs in
   our cohort for curated gene signatures for M1 or M2 macrophage states
   ([176]Fig. 4D). Consistent with similar single-cell analyses in other
   cancer types, GEP-NET TAMs did not separate neatly into M1 or M2
   groupings; for example, FOLR2-hi TAMs scored relatively highly for both
   the M1 and M2 signatures.

   We next sought to determine whether GEP-NET infiltrating TAMs displayed
   different chemokine profiles depending on the tumor site of origin. We
   performed a differential gene expression analysis within all TAMs in a
   given subtype grouped by tumor site of origin (i.e., pancreatic, small
   intestinal, or gastric) and examined the top five differentially
   expressed genes ([177]Fig. 4E). All TAMs originating from pNETs
   displayed increased expression of chemokines CCL3, CCL4, CCL4L2, and
   CCL3L1 ([178]Fig. 4, E and F). FOLR2-hi TAMs from siNETs differentially
   expressed chemokine CCL24, while FOLR2-hi TAMs and SELENOP-hi TAMs in
   the gNET sample had the highest expression of CXCL9 and CXCL10
   ([179]Fig. 4, E and F). CCL3, CCL4, CXCL9, and CXCL10 are known to
   recruit activated T cells via binding to the CXCR3 receptor, suggesting
   that FOLR2-hi and SELENOP-hi TAM subtypes identified in our cohort may
   have a role in orchestrating the antitumoral immune response in pNETs
   and gNETs. CCL24 has been implicated in activation of the M2 macrophage
   phenotype; SELENOP-hi TAMs in our siNETs may therefore contribute to a
   protumoral microenvironment. In summary, our results demonstrate the
   heterogeneity within the myeloid compartments of GEP-NETs and highlight
   promising therapeutic vulnerabilities common across GEP-NET subtypes.

Charting the transcriptional heterogeneity and evolution dynamics between a
primary siNET and its metastasis

   Our cohort contained a matched primary-metastasis pair from a patient
   with siNET, which provided a rare opportunity to investigate the
   temporal evolution of GEP-NETs. We separated cells from the primary
   tumor (sinet2) and its liver metastasis (sinet3) for combined analysis
   ([180]Fig. 5A). Both tumors contained a similar diversity of immune
   cells despite their different anatomical locations ([181]Fig. 5A). We
   then sought to compare the transcriptomic landscapes of the paired
   samples. Selection and unsupervised reclustering of only malignant
   cells from both tumors revealed four tumor clusters ([182]Fig. 5B).
   Sinet2, the primary tumor, comprised mostly cells from cluster 0 (C0),
   with a minority of cells distributed throughout the other three
   clusters ([183]Fig. 5C). Conversely, malignant cells from sinet3 fell
   predominantly in C1 and C2; no cells from the metastasis were found in
   C0 ([184]Fig. 5C).

Fig. 5. scRNA-seq reveals intratumoral heterogeneity and evolution dynamics
in GEP-NETs.

   [185]Fig. 5.
   [186]Open in a new tab

   (A) UMAP of malignant and immune cells from a primary siNET (sinet2)
   and its paired liver metastasis from the same patient (sinet3). Cells
   are colored by cancer or immune cell type. Tumor cells are represented
   by the clusters within the dashed oval. (B) UMAP of only tumor cells
   from sinet2 and sinet3, colored by Louvain cluster. (C) Donut plots
   representing percentage of cells found in each cluster from (B) for
   sinet2 and sinet3. Plots are colored by clusters as in (B). (D) Heatmap
   of scaled, normalized, and averaged expression of the top 10
   differentially expressed genes between tumor clusters from (B). (E)
   Copy number alteration profiles of sinet2 and sinet3. Tumor subclones
   are indicated by color bars on the left-hand side of each plot. (F)
   GSEA of epithelial-mesenchymal transition (top) and KRAS signaling up
   (bottom) transcriptional signatures in clone 1 cells compared to clone
   2 cells from sample sinet2.

   To interrogate the differences in transcriptional programming between
   each tumor cluster, we next conducted differential gene expression
   analysis on the malignant cells from sinet2 and sinet3 (Methods). We
   identified high expression of immediate early genes such as JUN, FOS,
   and EGR1 in C0; expression of these genes was also identified in tumor
   C3 ([187]Fig. 5D). C1 and C2 were distinguished from the other two
   clusters by differential expression of neuroendocrine genes such as
   CHGA and SLC18A1. Given that C1 and C2 contained mostly cells from
   sinet3, these data reflect the course of metastatic evolution of
   sinet2.

   To better understand the differential genetic alterations between the
   primary tumor and its metastasis, we inferred copy number variation
   profiles for both tumors. The primary tumor (sinet2) contained two
   clones, with clone 1 distinguished from clone 2 by amplifications in
   chromosomes 1 and 8 and deletions in chromosome 2 ([188]Fig. 5E). Next,
   we sought to interrogate differences in genetic programming between the
   subclones within sinet2 using gene set enrichment analysis (GSEA;
   Methods). Clone 1 demonstrated a high enrichment score for the HALLMARK
   “Epithelial-Mesenchymal Transition” geneset, implicating this subclone
   as a potential cellular origin for the liver metastasis ([189]Fig. 5F).
   Clone 2 was enriched for the HALLMARK “KRAS Signaling Up” geneset
   ([190]Fig. 5F). While KRAS mutations have been identified in high-grade
   pancreatic neuroendocrine carcinomas, further research is required to
   understand the potential role of KRAS signaling in GEP-NET oncogenesis
   and progression ([191]37). Together, our findings further support the
   observation of metastatic transcriptional reprogramming in NETs.

DISCUSSION

   To elucidate the fundamental cellular compositions and genetic
   programming of GEP-NETs, we interrogated single-cell transcriptomes
   from eight tumors spanning three anatomical sites of origin. Within our
   findings are reflections of known biology of NETs. For example, pNETs
   uniformly demonstrated high differential expression of endocrine
   pancreatic genes and displayed more intertumoral heterogeneity compared
   to other NET subtypes. siNETs expressed genes associated with monoamine
   synthesis and neuroendocrine vesicle biology, consistent with the known
   tendency of these tumors to secrete serotonin. Through scRNA-seq
   analysis, however, we were also able to further elucidate the
   transcriptional and GRN heterogeneity of GEP-NETs. Across our GEP-NET
   samples, we identified differentially expressed regulons with
   established roles in neuroendocrine cell differentiation and fate
   determination that expressed across various time points of
   gastrointestinal development. We further found that the transcriptional
   profiles of malignant pNET cells resemble more than one cell type in
   the mature human pancreatic islet, including somatostatin-producing
   delta cells and pancreatic polypeptide–producing gamma cells. These
   data imply the potential for lineage plasticity in GEP-NETs or an
   endocrine precursor cell provenance for these tumors, two possibilities
   that should be explored in future biological investigations.

   The evolution of NETs as they metastasize remains poorly understood.
   While our sample size was limited, we were able to identify two
   subclones within one of our primary siNETs, one of which bore
   transcriptomic and copy number alteration similarities to the tumor
   cells in the patient’s liver metastasis. Cells with robust expression
   of the neuroendocrine genes CHGA and SLC18A1 comprised less than 15% of
   tumor cells at the primary tumor site but 100% of cells in the
   metastasis, suggesting that a phenotypic shift occurred over the course
   of tumor progression and metastasis. Our study is the first to
   demonstrate the transcriptional dynamics of metastatic NETs at
   single-cell resolution; future studies are warranted to investigate the
   functional and translational significance of the transcriptional
   programs we identified in our work.

   Immunotherapy has revolutionized cancer treatment over the past few
   decades and has demonstrated benefit across several solid tumor types.
   Recent clinical trials investigating immunotherapy in neuroendocrine
   cancers have focused on targeting classical checkpoint genes such as
   CTLA4 and the PD1-PDL1/PDL2 axes in the lymphoid compartment ([192]38,
   [193]39). While many of these studies are ongoing, preliminary results
   suggest that undifferentiated, higher-grade neuroendocrine neoplasms
   may respond to checkpoint inhibition, but that use of checkpoint
   inhibitors in well-differentiated NETs, like those in our study, would
   yield limited therapeutic benefit ([194]40–[195]43). Consistent with
   these clinical observations, within all GEP-NET cancer cells in our
   dataset, we identified little to no expression of currently targetable
   immune checkpoint ligands, highlighting the limited treatment options
   in this disease. Similarly, within the lymphoid compartment of
   GEP-NETs, we further identified relatively low expression of common
   immunosuppressive ligands and receptors across T and NK cell types.
   4-1BB-hi CD8^+ T cells and FGFBP2^− NK cells from our gNET exhibited
   the highest expression of immune checkpoint genes; however, this sample
   was the only higher grade (WHO G3) tumor of our cohort and the only one
   to receive somatostatin analog therapy, everolimus, and chemotherapy,
   all factors that could have contributed to this finding. Further
   translational studies will be required to determine the clinical
   utility of targeting lymphoid cells with immunotherapy in
   well-differentiated NETs.

   In contrast, within the myeloid compartment, we observed high levels of
   immunosuppressive gene expression in multiple cell types across GEP-NET
   sites of origin and within both primary and metastatic tumors (fig.
   S5). Our data highlight opportunities for therapeutic strategies in
   these rare cancers based on targeting alternative pathways in the
   myeloid microenvironment of GEP-NETs. For instance, high expression of
   VSIR across myeloid cells in GEP-NETs suggests that targeting the VISTA
   protein could rearm T cells and macrophages and thereby reinvigorate
   antitumoral immunity in these patients. Preclinical studies of the
   VISTA pathway in multiple solid tumor models (e.g., pancreatic cancer,
   melanoma, colorectal cancer, renal cell carcinoma, and glioma) have
   affirmed its role in tumoral immune evasion ([196]33, [197]44,
   [198]45). Early clinical trial data are assessing the safety profiles
   of anti-VISTA immunotherapies in solid tumors that also target the
   PD-1/PD-L1 checkpoint; per our findings, well-differentiated NETs may
   benefit from this therapeutic approach and should be considered for
   inclusion in future trial cohorts ([199]46, [200]47). Similarly, we
   identified expression of HAVCR2 (TIM3) and LGALS9 (Gal-9) in myeloid
   cells of the GEP-NET microenvironment. Combination anti-Tim3/anti–PD-1
   or anti–Gal-9/anti–PD-1 therapies have proven efficacious in enhancing
   T cell activation in preclinical studies of solid tumors and are
   currently under investigation in clinical trials ([201]48–[202]50).
   Last, inhibitors of SIGLEC10-CD24 binding between TAMs and cancer
   cells, respectively, are an emerging focus of preclinical investigation
   ([203]51). Our study therefore provides evidence supporting the use of
   these emerging immunotherapies, which may one day augment the limited
   effective treatment options available to patients with nonresectable
   GEP-NETs.

   Our study has several limitations. While the diversity of our cohort
   permits us to assess shared transcriptional programs across GEP-NET
   subtypes, our sample size remains relatively small due to the rarity of
   these tumors. Treatment history at the time of surgical resection also
   varied between cases; all pNETs and one of the siNETs were treatment
   naïve, while the metastatic small intestinal and gNETs had received at
   least somatostatin analog therapy. Differences in cellular phenotypes
   and microenvironment composition may therefore reflect therapeutic
   pressures and not necessarily innate biology. To assess these
   possibilities, future studies should interrogate large, clinically
   diverse tumor cohorts that include GEP-NET subtypes and treatment
   modalities not represented in our study. Inclusion of more
   treatment-naïve cases can also help to better define innate tumor
   biology. Last, we here present an analysis of a singular
   well-differentiated, G3 gNET, a classification defined by the highest
   mitotic index out of all well-differentiated NETs. Given the high
   degree of similarity between well-differentiated G3 GEP-NETs and poorly
   differentiated, aggressive neuroendocrine carcinomas, identification,
   and inclusion of G3 GEP-NETs in studies of well-differentiated tumors
   demands careful review by experienced pathology teams. Despite the
   heterogeneity of cases included in our study, our cohort encompasses
   multiple, well-differentiated GEP-NET subtypes across diverse treatment
   histories; our findings therefore serve as a key foundation for future
   basic and clinical investigations into these rare, heterogeneous
   tumors. We anticipate that in addition to further evaluating our
   finding of common immunosuppressive myeloid phenotypes, new axes of
   shared biology can be identified through future multimodal sequencing
   of the transcriptional and epigenomic landscapes of GEP-NETs.

   From a technical perspective, although we were able to infer copy
   number alteration profiles for each tumor using the inferCNV approach,
   verification of these inferred profiles will require patient-matched
   whole-genome sequencing. Overall, inferring tumor subclone architecture
   from single-cell data remains a challenging problem for the field given
   the sparsity and sampling biases inherent to scRNA-seq. Future work
   that obtains matched scRNA-seq and whole-genome sequencing can robustly
   characterize malignant subclones and also examine the influence of a
   GEP-NET patient’s genetic background on their tumor’s cellular
   landscape.

   In summary, our findings uncover previously unappreciated heterogeneity
   within NET subtypes and reveal potential evolution in tumor
   characteristics as they metastasize. Consistent with other studies, we
   also found that in well-differentiated NETs, expression of checkpoint
   markers targetable with current standard immunotherapies is uncommon.
   However, we also identified immunosuppressive markers in myeloid cells
   that are commonly expressed in NETs and may form the basis for future
   targeted therapies. Our findings therefore provide an important step in
   delineating key underlying molecular features of NETs and shed light on
   potential novel therapeutic strategies for patients suffering from this
   challenging disease.

METHODS

Experimental model and subject details

   Human patient samples were collected under Dana-Farber Cancer Institute
   IRB protocol #02-314. Informed consent was required for this study and
   was obtained as follows: At the time of their clinic appointments,
   patients were contacted by a member of the study staff and asked if
   they were interested in participating in this research study involving
   collection of biological samples and clinical data in patients with
   NETs. The study was explained in detail to interested patients, and
   his/her questions were answered. Informed consent was signed to
   participate. Patients were specifically consented to (i) allow the
   collection of tissue, blood, or urine; (ii) to allow the collection and
   linkage of clinical data to the tissue; (iii) to allow for recontact of
   the patient for the purpose of recruitment into new clinical protocols,
   or for the purpose of obtaining clinical follow-up; (iv) to allow for a
   one-time additional blood draw or urine collection for study purposes
   if not feasible in the context of a scheduled blood draw; and (v)
   collection of periodic blood samples. The consent status of each
   patient was recorded in the database, which allows specific tailoring
   of data acquisition according to each patient’s wishes. These measures
   will ensure that the wishes of each patient are maintained.

Methods details

Sample collection and dissociation for scRNA-seq

   Fresh tissue from surgical resections was collected at Brigham and
   Women’s Hospital. Members of the pathology team evaluated surgical
   specimens to identify tumor tissue. Tumor tissue that was not required
   for routine clinical pathological assessment was dissected away from
   nontumor tissue using a scalpel and transferred to a 50-ml conical tube
   with Hank’s balanced salt solution and transported on ice. On arrival
   to the laboratory, tumor samples were transferred into a 5-ml Eppendorf
   tube containing 4.5-ml cold enzymatic dissociation mix [collagenase
   type 4 (100 μg/ml; Worthington Biochemical Corporation, #[204]LS004186)
   and deoxyribonuclease I (100 μg/ml; StemCell Technologies, #07900) in
   RPMI + Hepes (Thermo Fisher Scientific, #22400089)]. The resection
   tissue was then minced inside the tube with spring scissors into
   <0.5-mm fragments at room temperature. The resulting mixture was
   incubated for 10 min in a 37°C water bath with manual inversions every
   1 to 2 min to ensure thorough mixing. Next, the mixture was vigorously
   and repeatedly pipetted using a P1000 pipette at room temperature for
   further mechanical dissociation. For all samples except pnet2 and
   pnet3, the tube containing the mixture was then returned to the 37°C
   water bath for an additional 10 min with manual inversions every 1 to 2
   min for thorough mixing.

   The cell suspension was then filtered through a 70-μm cell strainer
   into a 15-ml conical tube on ice. The strainer was washed with 5-ml
   cold RPMI to maximize cellular yield. The filtered cells then underwent
   centrifugation at 450g for 5 min at 4°C. The resulting supernatant was
   carefully transferred to a 15-ml conical tube on ice, taking care not
   to disturb the pellet. To lyse red blood cells, the centrifuged cell
   pellet was resuspended in 300 to 500 μl of ACK lysing buffer (Thermo
   Fisher Scientific, #A1049201) and incubated for 1 min on ice. Calcium-
   and magnesium-free phosphate-buffered saline (PBS; Thermo Fisher
   Scientific, #10010-23) was added in a volume twice that of the ACK
   lysing buffer to stop red blood cell lysis. Cells were transferred to a
   1.7-ml Eppendorf tube on ice and then recentrifuged for 8 s at 4°C with
   centrifugal force ramping up to but not exceeding 11,000g. After
   recentrifugation, cells were resuspended in calcium- and magnesium-free
   PBS with 0.4% bovine serum albumin (BSA; Ambion, #AM2616). If a large
   amount of red blood cells remained after recentrifugation, an
   additional round of red blood cell lysis was performed. After
   satisfactory red blood cell lysis, the supernatant was carefully
   transferred to a new 15-ml conical tube on ice, and the remaining cell
   pellet was resuspended in 50 μl of RPMI + Hepes.

   Five microliters of the single-cell suspension was mixed with Trypan
   blue (Sigma-Aldrich, #T8154) at a 1:1 volume ratio and loaded into a
   hemocytometer for cell counting and assessment of cell viability,
   clumping, and debris. At the end of processing, all samples had
   viability greater than 51%. Cell suspensions were then centrifuged and
   resuspended in either RPMI + Hepes or RPMI + Hepes with 0.4% BSA in
   preparation for transport and downstream flow cytometric processing.

Flow cytometry

   Following viability assessment, samples were centrifuged at 600g for 5
   min and resuspended in fluorescence-activated cell sorting (FACS)
   buffer [1% BSA (Millipore Sigma, product no. EM-2930) in
   calcium/magnesium-free PBS (Thermo Fisher Scientific, product no.
   14190144)] and stained with 1 μg of CD45-APC antibody (Thermo Fisher
   Scientific, product no. 50-166-056) for 30 min at 4°C. The samples were
   washed with 1 ml of FACS buffer and centrifuged at 600g for 5 min
   before being resuspended in 100 μl of FACS buffer with 10 μM
   Stemolecule Y27632 (ROCK inhibitor) (Stemgent, product no. 04-0012-02).
   A total of 200 μl of FACS buffer with calcein blue AM (Thermo Fisher
   Scientific, product no. C1429) was added to the sample for live/dead
   exclusion and the sample was passed through a 35-μm strainer before
   sorting. Two fractions of the samples (live and CD45^+) were sorted
   into FACS buffer using a Beckman Coulter MoFlo Astrios cytometer.

Single-cell RNA sequencing

   For each sample, the live bulk and CD45^+ fractions were centrifuged at
   600g for 5 min and resuspended in FACS buffer for a final concentration
   of 600 cells/μl as specified in the 10x Genomics user guide document
   number [205]CG000185, Rev. B for an estimated recovery of 3000 cells.
   These fractions were loaded into the 10x Genomics Single Cell chip (10x
   Genomics, #PN-1000120) along with gel beads and reverse transcription
   master mix from the Single Cell 5′ v1.1 Chemistry Kit (10x Genomics,
   #PN-1000165). The chip was run on the 10x Chromium Controller using the
   standard program to produce complementary DNA, which was then amplified
   and used for Gene Expression (GEX) library generation. All single-cell
   GEX libraries were sequenced using an Illumina NextSeq 2000 sequencer.

scRNA-seq data preprocessing

   Raw sequencing data underwent demultiplexing, barcode processing, read
   alignment, and UMI counting using the 10x Genomics CellRanger pipeline
   (version 5.0) with default parameters. Reads were aligned to the
   prebuilt human genome reference included in Cell Ranger (GRCh38). After
   filtering to exclude low-RNA content droplets, a gene-barcode matrix
   was generated for each sample containing counts of confidently mapped,
   non–PCR-duplicate reads.

   To remove technical artifacts arising from cell-free RNA profiles from
   our sequencing data, ambient RNA decontamination was performed on all
   samples using the CellBender software package ([206]52). In brief, the
   raw counts matrix and expected cell count for each sample were provided
   as input to the remove-background module of CellBender, which uses an
   unsupervised deep-generative model to distinguish empty droplets from
   cell-containing ones. Expected cell count was derived from the estimate
   generated by Cell Ranger. All samples were run using the “ambient”
   model and default settings of 150 training epochs and false-positive
   rate of 0.01. Transcripts assigned to empty cells were iteratively
   detected and removed from the raw counts matrix by remove-background,
   yielding a “background-subtracted” cleaned counts matrix as output.

   To remove multiplet artifacts caused by droplet encapsulation of more
   than one cell, we then performed doublet detection and removal on the
   CellBender-cleaned counts matrix using the Scrublet package (Python,
   version 0.2.2) ([207]53). An expected doublet rate of 0.06 and manually
   selected doublet score thresholds were used to separate singlets from
   neotypic doublets in the resulting bimodal score distribution graphs.
   Presumed doublets were removed to create an ambient RNA-free,
   singlet-only cleaned counts matrix for downstream analysis. All further
   quality control, dimensionality reduction, unsupervised clustering, and
   differential expression analyses were performed using the Seurat R
   package (version 4.0.6) ([208]54).

Cell type clustering and across-sample integration

   Before clustering, the 16 datasets (one unsorted and one CD45^+ sorted
   per sample) were merged into a single Seurat object. Cells with fewer
   than 200 genes or greater than 25% counts representing mitochondrial
   genes were removed from each sample. Genes detected in fewer than three
   cells per sample were also removed. The raw counts data were then
   normalized and scaled using the default log-normalization approach from
   the R package Seurat (version 4.0.6) ([209]54).

   Linear dimensional reduction was performed on the integrated and
   normalized counts matrix by running principal components analysis (PCA)
   on the 5000 most highly variable genes in our dataset. The first 30
   principal components were used to perform Louvain clustering with a
   resolution parameter of 0.4. Uniform Manifold Approximation and
   Projection (UMAP) was then performed using the same principal
   components to visualize preliminary cell clusters from the
   nonintegrated dataset in two-dimensional space.

   To mitigate technical and biological variation between samples, the
   merged dataset was integrated using the R package Harmony (version
   0.1.0), with patient of origin as the primary source of variation
   ([210]55). Harmony uses the PCA embeddings for each cell to iteratively
   generate multisample clusters, penalizing clustering arrangements that
   contain clusters with one or few samples represented. The final
   integrated dataset included 24,048 cells from the eight-tumor cohort.

Cell type visualization, and identification

   To determine cell types within each cluster, differential gene
   expression analysis was performed by comparing cells from each cluster
   to all other cells using a two-sided Wilcoxon rank-sum test with
   applied Bonferroni correction for multiple comparison testing. Top
   differential genes within each cluster were used to assign cells into
   broad tumor and immune categories. For subclassification of immune cell
   types, lymphoid and myeloid cells were subset from the integrated
   dataset and clustered and UMAP projected using 30 principal components.
   The same differential gene expression analysis was used to identify T
   cell and myeloid types at higher resolution. CD8^+ T cells and TAMs
   were further reclustered and reprojected to conduct similar cellular
   dissection.

Identification of cancer cells

   To identify malignant cells, we first integrated both non–CD45^+-sorted
   and CD45^+-sorted datasets from the same tumor to generate a combined
   expression matrix. Each sample underwent Louvain clustering and UMAP
   projection to identify PTPRC^+ (CD45^+) and PTPRC^− clusters. We then
   used inferCNV (R, Terra implementation) to estimate the CNV profile of
   each log-normalized sample using all PTPRC^+ cells as the reference
   group and all PTPRC^− cells as the observation group. NET-related
   marker gene expression and presence of inferred copy number alteration
   as determined by inferCNV were combined to identify “tumor” cells
   within our cohort.

Differential expression, GSEA, and signature scoring

   To minimize the effect of patient-specific expression patterns on our
   differential gene expression analysis of GEP-NET cells across
   anatomical subtypes ([211]Fig. 2C), we averaged the expression of genes
   from all tumors from the same primary site of origin (i.e., pancreas or
   small intestine) and compared the two transcriptional profiles. To
   compare the transcriptomic profiles of cells from a paired primary
   (sinet2) and metastatic (sinet3) tumor set, differential gene
   expression analysis was performed using a two-sided Wilcoxon rank-sum
   test with Bonferroni correction applied. The log[2](fold change) values
   for each gene were used as ranks for preranked GSEA, which was
   performed using the fgsea R package (version 1.20.0) ([212]56). GSEA
   was performed using all HALLMARK, Kyoto Encyclopedia of Genes and
   Genomes, and Gene Ontology Biological Process gene sets from version
   6.2 of the MSigDb as input ([213]57). Significance level was set at P <
   0.05 and false discovery rate < 0.25.

   Single-cell signature scoring was performed using the VISION R package
   (version 3.0.0) ([214]58). Progenitor exhaustion and terminal
   exhaustion signatures derived from Bi et al. ([215]30) were used to
   score CD8^+ T cells from all tumors ([216]Fig. 3F). Known macrophage,
   monocyte, and DC markers were used to generate custom gene signatures,
   which were then used for VISION scoring of the NET myeloid compartment
   ([217]Fig. 3B).

SCENIC analysis and regulon inference

   To identify differentially active GRNs in tumor cells across GEP-NETs,
   we performed SCENIC analysis using the recommended pipeline provided
   with the pySCENIC implementation
   ([218]https://github.com/aertslab/SCENICprotocol/blob/master/notebooks/
   PBMC10k_downstream-analysis.html) ([219]59, [220]60). We first
   performed GRN inference using the multiprocessing implementation of
   Arboreto to run GRNBoost2 on the count matrix containing averaged
   expression of tumor cells from each sample. We then ran the SCENIC
   implementation of cisTarget to predict regulons from the TF
   adjancencies matrix produced by GRNBoost2. Last, we used AUCell to
   generate regulon specificity scores for each tumor. As recommended in
   the pySCENIC documentation, we selected the top five regulons in each
   tumor ranked by AUC, then plotted the scaled scores for these regulons
   by tumor.

Metamodule scoring

   To generate metamodule scores for TAMs in [221]Fig. 4, we first used
   Seurat’s FindAllMarkers function to identify TAM cluster–specific
   marker genes. We then used the average log[2](fold change) of genes in
   our differential expression analysis as a ranking metric for GSEA
   analysis on all TAMs using the HALLMARK genesets from the MSigDB. Four
   of the most enriched pathways from this analysis were selected as
   biologically relevant metamodules. To score each TAM for expression of
   these four metamodules, we implemented the method developed in Neftel
   et al. ([222]34) using the scrabble R package. In brief, a signature
   score (SC) is calculated for each cell by averaging the relative
   expression of the signature genes in that cell, then subtracting the
   average relative expression of a control geneset in the same cell. The
   control geneset is defined as previously ([223]34). All cells were
   scored for each of the four metamodules. The position of each cell in
   the butterfly plot ([224]Fig. 4C) was determined by calculating
   [sign(SC1 − SC2)
   [MATH: <mo>∗</mo> :MATH]
   log[2](|SC1 − SC2| + 1)] for both x and y coordinates.

Determination of CNV clonal structure

   To determine the clonal variation between the primary (sinet2) and
   metastatic tumor cells (sinet3) from the same patient with a siNET,
   cells from both samples were processed using the pipelineCNA function
   from the R package SCEVAN ([225]61). The “subclones” parameter was set
   to “TRUE” to enable analysis of clonal structure within each tumor.
   Heatmaps displaying inferred copy number loss/gain across chromosomes 1
   to 22 and showcasing the number of inferred subclones were generated
   for each sample through the SCEVAN pipeline. To identify pathways
   enriched in clone 1 and clone 2 from the primary tumor (sinet2),
   differential gene expression was performed comparing the pseudobulked
   malignant cells from each clone using Seurat’s FindAllMarkers function.
   The resulting differential gene list for each clone was ranked by log2
   fold change and used as input for GSEA using the fgsea package (R,
   version 3.17).

Multiplex immunofluorescence staining and imaging

   Multiplex immunofluorescence staining was performed on five FFPE
   samples from our cohort using a BOND RX fully automated stainer (Leica
   Biosystems). FFPE tissue sections of 5 μm thick are baked for 3 hours
   at 60°C before loading into the BOND RX. Slides are deparaffinized
   (BOND DeWax Solution, Leica Biosystems, catalog no. AR9590) and
   rehydrated with series of graded ethanol to deionized water. Antigen
   retrieval is performed in BOND epitope retrieval solution 1 (pH 6) or 2
   (pH 9), as shown in table S2 (ER1 and ER2, Leica Biosystems, catalog
   nos. AR9961 and AR9640) at 95°C. Deparaffinization, rehydration, and
   antigen retrieval are preprogrammed and executed by the BOND RX. Slides
   are then serially stained with primary antibodies. Incubation time per
   primary antibody was 30 min. Subsequently, anti-mouse plus anti-rabbit
   Opal Polymer Horseradish Peroxidase (Opal Polymer HRP Ms. + Rb, Akoya
   Biosciences, catalog no. ARH1001EA) is applied as a secondary label
   with an incubation time of 10 min. Signal for antibody complexes is
   labeled and visualized by their corresponding Opal Fluorophore Reagents
   (Akoya Biosciences) by incubating the slides for 10 min. Opal
   Fluorophore 780 is paired with a TSA-DIG amplification to ensure
   adequate and analyzable signal. Slides are incubated in Spectral DAPI
   solution (Akoya Biosciences) for 10 min, air dried, and mounted with
   Prolong Diamond Anti-fade mounting medium (Life Technologies, catalog
   no. [226]P36965) and stored in a light-proof box at 4°C before imaging.
   The target antigens, antibody clones, dilutions for markers, diluents,
   and antigen retrieval details are listed in table S2.

   Image acquisition is performed using the PhenoImager multispectral
   imaging platform (Akoya Biosciences, Marlborough, MA). Each slide is
   scanned at 20× resolution as whole-slide scan images. These images are
   then accessible through Phenochart viewing software (Akoya Biosciences)
   where four to eight 20× regions of interest (ROIs) are selected. On the
   basis of well-staining tissue availability, three to seven ROIs per
   sample were chosen for analysis by a board-certified anatomic and
   molecular pathologist specialized in tissue-based biomarker assays.
   ROIs were chosen based on high tumor content; areas of necrosis,
   fibrosis, and unusual vascular architecture were explicitly avoided.

   After ROI selection and approval from a pathologist (S.R.), the images
   are spectrally unmixed and analyzed within Inform 2.6 (Akoya
   Biosciences). The ROIs are subsequently segmented and quantified for
   expression of each marker using the Inform analysis tools. Data tables
   are exported from Inform and run through a custom data extraction
   pipeline to obtain cell population densities (number of cells per
   square millimeter) for each marker and/or combinations of markers.

Acknowledgments