Gene Expression Signatures Identify Novel Therapeutics for Metastatic Pancreatic Neuroendocrine Tumors

Pancreatic neuroendocrine tumors (pNETs) are uncommon neoplasms that account for less than 3% of all pancreatic tumors (1). However, the annual incidence of all neuroendocrine tumors (NET) has increased steadily over the past several decades, and pNETs now have an approximate annual incidence of 0.8 per 100,000 persons (2). They are characteristically slow-growing malignancies associated with a comparably favorable prognosis, especially when contrasted against the significantly more common pancreatic adenocarcinoma. The indolent nature of the disease often delays diagnosis such that metastases, predominantly to the liver, are present in approximately 60% of patients at diagnosis (3, 4). Despite advances in knowledge and management of NETs, the majority of patients with metastatic pNETs will die of their disease, and current therapies have not been shown to definitively improve overall survival. Greater understanding of the biology of these tumors, particularly the changes associated with metastasis, is needed to identify new therapeutic targets and more effective treatments.

Clinical management of metastatic pNETs is multimodal and rapidly evolving. Historically, surgical excision of the primary tumor, debulking of metastatic lesions, and somatostatin analogues (SSA) were the main therapeutic options. While surgery remains the standard of care for localized tumors and an important component of the treatment of metastatic disease, several systemic agents have emerged for treating pNETs. Streptozocin (STZ) was one of the first drugs found to be active against pNETs, although low response rates and significant toxicity limit its use today (5). SSAs have long been used to reduce the symptoms of excess hormone secretion associated with many NETs, and have more recently been shown to have antiproliferative effects as well (6, 7). Currently, SSAs are considered the first-line therapy for metastatic NETs expressing somatostatin receptors. In the past decade, the options for treatment have expanded further with evidence of improved progression-free survival (PFS) using the tyrosine kinase inhibitor, sunitinib, and the mTOR inhibitor, everolimus (8, 9). Other tyrosine kinase inhibitors, including cabozantinib and sulfatinib, have shown promise in phase II trials, and phase III trials are underway (10). The combination of the antimetabolite capecitabine and alkylating agent temozolomide (CAPTEM) can be given orally and has less toxicity as compared with STZ-based treatment (11–13). A phase II study has shown significant PFS benefits of this combination over temozolomide alone (14). Peptide receptor radionuclide therapy (PRRT), which uses a radiolabeled SSA to target NETs, was recently approved for the treatment of gastroenteropancreatic NETs after impressive improvements in PFS of treated patients with grade 1 and 2 small bowel NETs relative to those receiving a high-dose SSA (15). Therapies that extend overall patient survival, however, are still needed.

Pancreatic neuroendocrine tumors have been recognized as a distinct clinical entity for nearly a century, yet our understanding of their genetic basis is incomplete. Several hereditary cancer syndromes are associated with pNET development, including multiple endocrine neoplasia type 1 (MEN1), von Hippel-Lindau disease (VHL), neurofibromatosis type 1 (NF-1), and tuberous sclerosis complex (TSC1 and 2), but the majority of pNETs occur sporadically (>90%; ref. 16). Among sporadic pNETs, the most commonly altered gene is MEN1, with mutations found in 41%–44% of cases. Other commonly mutated genes include DAXX or ATRX, genes in the mTOR pathway, genes involved in DNA damage repair such as MUTHY, and genes involved in chromatin modification (17, 18). Although mutations in the mTOR pathway are seen in only 12%–15% of pNETs, expression and activity level analyses have revealed alterations of PI3K/Akt/mTOR signaling involving both upstream and downstream regulators in the majority of these tumors (16, 19, 20).

In this study, we investigated gene expression patterns in nodal and liver metastases versus primary pNETs (most of which were patient-matched) to identify drug-targetable regulators of metastasis. We uncovered prominent alterations in MAPK, cyclin-dependent kinase (CDK), topoisomerase (TOP2A), and NF-kB signaling along with expected changes in PI3K and mTOR pathways. Analyses of the metastatic tumor transcriptional changes using two complementary bioinformatic tools identified drugs that displayed significant in vitro efficacy against pNET cell proliferation and survival. The approach, identification of NET metastatic pathways, and discovery of new anti-NET drugs described herein offers unique insights and promising therapies for treating advanced pNETs.

Tumor samples were obtained from 43 patients undergoing resection of well-differentiated, grade 1 and 2 pNETs between December 11, 2006 and August 19, 2016 (Supplementary Table S1). Tissue samples from the primary tumor, liver, and lymph node metastases were collected from all patients where available. Informed written consent was provided by all patients in accordance with a protocol approved by the University of Iowa Institutional Review Board (IRB Number 199911057), and studies were conducted in accordance with the Belmont Report. Following surgery, tissue samples were stored at −20°C in RNALater (Thermo Fisher Scientific, Waltham). Total RNA was extracted using the RNeasy Plus Universal Mini Kit (Qiagen) per the protocol recommended by the manufacturer. Tumor cellularity was estimated by a pathologist with NET expertise (A.M. Bellizi) from paraffin-embedded tissue stained with hematoxylin and eosin, and samples with less than 70% neuroendocrine cellularity were excluded. The Agilent 2100 Bioanalyzer (Agilent Technologies) was used to assess RNA quality and assign RNA integrity numbers (RINs).

Transcription profiling using RNA-Seq was performed by the University of Iowa Genomics Division using manufacturer-recommended protocols. Initially, 500 ng of DNase I-treated tRNA was used to enrich for polyA-containing transcripts using oligo(dT) primers bound to beads. The enriched RNA pool was then fragmented, converted to cDNA, and ligated to sequencing adaptors containing indexes using the Illumina TruSeq stranded mRNA sample preparation kit (catalog No. RS-122-2101, Illumina, Inc.). The molar concentrations of the indexed libraries were measured using the 2100 Agilent Bioanalyzer (Agilent Technologies) and combined equally into pools for sequencing. The concentrations of the pools were measured using the Illumina Library Quantification Kit (KAPA Biosystems) and sequenced on the Illumina HiSeq 4000 genome sequencer using 150-bp paired-end sequencing by synthesis (SBS) chemistry.

FASTQ files were used to quantify transcript abundance using Salmon version 0.11.3 (21). Tximport 1.8.0 was used to perform gene-based quantification of abundance, using Ensembl Genes version 96/Biomart to map gene symbols to Ensembl transcript identifiers (22, 23), and voom (24) for library normalization and to generate normalized transcripts per million (TPM) values. StringTie version 1.3.0 (25) was used to assign reads to transcripts using the GFF3 file from Ensembl (http://ftp.ensembl.org/pub/grch37/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh37.87.gtf.gz).

RNA-seq count data was normalized using TPM. To identify differentially expressed genes, the nonparametric Wilcoxon rank-sum test was used. Log₂-fold changes (LFC) between tumors at different sites were estimated using the voom package (24). LFC > 0 indicates that a gene was upregulated in the metastases, whereas LFC < 0 indicates downregulation. To filter out results that may result from the confounding influence of comparing different tissue types between primary and metastatic tumors, we carried out two rank-sum tests: one comparing primary and metastatic tumors (combining liver and lymph node metastases), and the other comparing liver and lymph node metastases. From the results of these two tests, genes of interest are those that are significantly differentially expressed between primary and metastatic sites (P < 0.05), but not differentially expressed between liver and lymph node tumors. Genes were defined as tissue-specific and filtered out if expression was twofold higher at either metastatic site compared with the other (LFC > 1). This tissue-specific filtering was designed to remove genes whose overexpression was dependent on the site of metastasis, and thus was applied only to the genes that were upregulated in the metastases. All analyses were performed in R 3.5.0 (www.R-project.org). Unsupervised analysis of the PNET samples was performed using both t-distributed stochastic neighbor embedding (t-SNE) and principal component analysis (PCA).

The filtered list of differentially expressed metastatic genes was used to query the Connectivity Map touchstone database (CMap; Broad Institute, ref. 26). This database includes gene expression signatures derived from nine cancer cell lines treated with 2,429 well-annotated compounds. These expression profiles were compared with our metastatic PNET signature, and compounds were assigned a connectivity score (ranging from −100 to 100) based on the similarity of the expression changes they induced when compared with our metastatic signature. A positive connectivity score indicates similar gene expression changes, while a negative score indicates an opposing pattern. Each compound and cell line are assigned a separate score, and a summary score is provided for each compound across all cell lines. In parallel, the IPA analysis match function was also used to predict potentially therapeutic drugs. The IPA output of our metastatic profile was compared with analyses from other cell lines treated with various compounds, and these drugs were assigned a score that ranged from −14.3 to −244.3. A more negative score indicated expression changes resulting in predictions that were in the opposite direction of those based on our metastatic profile. The molecule activity predictor and upstream analysis tools in IPA (Qiagen; ref. 27) were used to predict the activation state of drug targets from our metastatic PNET expression signature. IPA calculates a z-score to predict the activation state of upstream regulators, with a positive score indicating activation, a negative score indicating inhibition, and an absolute value greater than 2 suggesting significance. Drugs were selected for testing in pNET cell lines based on a combination of factors: a highly negative connectivity score from CMap and the IPA analysis match, drug availability, and a mechanism of action that was highly represented on both the CMap and IPA lists.

A subset of 14 significantly differentially expressed genes were selected for qPCR validation in a separate cohort of PNET primary and metastatic samples not tested by RNA-Seq. RNA was extracted as described above from primary tumors, lymph node metastases, and liver metastases from 12 other available patients with well-differentiated, grade 1 or 2 pNETs. Total RNA was then reverse-transcribed to cDNA using the qScript cDNA Supermix (QuantaBio). qPCR was performed using gene-specific primers from Integrated DNA Technologies and PerfeCTa SYBR Green Supermix Dye (QuantaBio) using the 7900HT Fast Real-Time PCR System (Applied Biosystems). Primer sequences were obtained from PrimerBank (https://pga.mgh.harvard.edu/primerbank/). Genes were selected using IPA from a list of genes whose differential expression was used to predict the activation states of our drug targets using the molecule activity predictor. Expression was normalized against the internal control gene RPLP1, which was used to calculate the delta cycle threshold (dCT) for each gene of interest. The average expression of each gene in the metastatic (liver and lymph node) tumors was compared with expression in the primary tumors, and statistically significant differential expression was determined using the Wilcoxon rank-sum test.

Drugs were purchased from Medchem Express LLC and Thermo Fisher Scientific, dissolved in DMSO, and stock solutions stored at −20°C. BON-1 and QGP-1 pNET cells were seeded at a density of 2,000 cells per well in 96-well flat-bottomed culture dishes. After overnight incubation, each drug was added at the indicated concentrations and incubated for 5 days; assays were performed in triplicate. Samples were evaluated for relative cell number using AlamarBlue (Thermo Fisher Scientific). Results were quantified using a fluorescence microplate reader by measuring fluorescence of AlamarBlue at an excitation wavelength of 560 nm with fluorescence emission at 590 nm. For simvastatin-treated cells, media containing simvastatin was removed and replaced with new media prior to the addition of the AlamarBlue reagent because high concentrations of simvastatin interferes with measurement of AlamarBlue. Results were analyzed using CompuSyn Software (ComboSyn, Inc.) to determine the IC₅₀ for each drug.

Two human PNET cell lines were used for these studies. BON-1 cells, established from a lymph node metastasis of a pNET (28), were maintained in DMEM/F12 containing 10% FBS, 4 mmol/L glutamine, 100 U/mL penicillin, and 100 μg/mL streptomycin. QGP-1 cells, established from primary somatostatin producing pNET (29), were purchased from the Japanese Collection of Research Bioresources (JCRB0183), and were maintained in RPMI1640 medium containing 10% FBS, 4 mmol/L glutamine, 100 U/mL penicillin, and 100 μg/mL streptomycin. Both cell lines were recently authenticated by Hofving and colleagues and found to express neuroendocrine markers and harbor mutations associated with NETs (30). Both cell lines were tested for Mycoplasma contamination by selective enzymatic assay (MycoAlert Plus, Lonza Inc.) and found to be negative (most recently on January 6, 2020). Cryopreserved stocks of each cell line were made at low passage numbers (less than 5–8 passages after receipt) and cells were used in experiments for a maximum of 10–12 passages (roughly 5–6 weeks) after being thawed.

RNA-Seq was performed on 39 primary tumors, 21 lymph node metastases, and 17 liver metastases from 43 patients with pNETs (Supplementary Table S1). All patients had grade 1 or 2, well-differentiated tumors. A flowchart summarizing our methods from sequencing to drug testing is presented in Fig. 1. A total of 18,103 genes were measurably expressed in the tumor samples. The results of our RNA-Seq are available on the database of Genotypes and Phenotypes (dbGaP) under phs study number phs001772.v1.p1 (https://www.ncbi.nlm.nih.gov/projects/gapprev/gap/cgi-bin/preview1.cgi?GAP_phs_code=DrKURjFTCAuySlUS). Of these, 902 genes were significantly differentially expressed in the metastatic tumors compared with the primaries: 718 were overexpressed and 184 underexpressed (Fig. 2A). Of the overexpressed genes, 13 showed greater than twofold higher expression in the nodal versus liver metastases, while 263 were expressed at greater than twofold higher levels in the liver versus nodal metastases. These 276 genes were filtered out to reduce tissue-specific effects on the metastatic gene signature, leaving 626 differentially expressed genes, which constituted our final metastatic pNET gene signature (Supplementary Table S2; Fig. 2B).

Analysis of the liver genes that were removed by the filter using IPA revealed that they were largely involved in processes associated with normal liver function, including production of coagulation factors and lipid metabolism. Moreover, IPA analysis of the unfiltered metastatic gene set (902 genes) demonstrated the dramatic effect that inclusion of these liver genes had on subsequent analysis; in the unfiltered set, the pathways and cellular processes with the strongest predicted activation were generally associated with normal hepatic function.

Dimensionality reduction using PCA and t-SNE revealed that the primary determinant of clustering was the tissue of origin (Fig. 3). This was particularly true for the liver metastases that tended to cluster separately from the primary tumors and lymph node metastases, which were much closer together. Despite this tendency, there were several outliers for which the clustering was primarily determined by the individual patient. This second trend is most clearly visualized by t-SNE, which also shows that six of the liver metastases cluster tightly with their corresponding primary and lymph node samples, rather than with the other liver metastases.

We queried the CMap touchstone database with expression signatures from the filtered metastatic PNET gene expression profile. The top 100 compounds producing expression changes that were the reverse of those seen in our metastatic profile are shown in Supplementary Table S3. Among the drugs with the lowest connectivity scores, inhibitors of mTOR and PI3K were the most highly represented. In parallel, the analysis match functionality of IPA was used to identify compounds that featured related pathways with an opposing activation state for prioritization. The most highly rated compounds predicted by IPA analysis match are shown in Supplementary Table S4, with numerous inhibitors of mTOR, PI3K, BRAF/MEK, and CDK ranked highly. Fifteen compounds were selected from these two lists using the criteria described previously (Table 1). These compounds included inhibitors of CDK, histone deacetylase (HDAC), NF-kB signaling, DNA topoisomerase, RNA-polymerase, DNA protein kinase, MEK, hydroxymethylglutaryl-coenzyme A reductase (HMGCR), and the PI3K/AKT/mTOR pathway.

The upstream analysis tool in IPA was used to analyze the metastatic gene expression signature and a number of regulators were predicted to be activated using a z-score cutoff of 2: OSCAR, MEK, IL27, NFE2L2, CD40LG, CAMP, TREM1, PRKCD, collagen type II, and NANOG (Supplementary Table S5). Three upstream regulators were also predicted to be inhibited with z-scores less than −2: CST5, miR-146, and PDCD1. With the exceptions of MEK and protein kinase C delta (PRKCD), none of the other drug targets suggested by CMap were independently predicted to be activated by the IPA upstream analysis tool. In comparison, the molecule activity predictor tool in IPA was used to make predictions as to the activation state of several targets of the drugs predicted by CMap, using a more expansive set of connections that includes data derived from nonhuman studies. Using this analytic approach, CDK, NF-kB, mTOR, MEK, HDAC1, PI3K, and protein kinase C (PKC) were all predicted to be activated in the metastases. These predictions are shown in graphical form in Fig. 4.

Fifteen genes found to be differentially expressed in our RNA-Seq data were selected for validation using qRT-PCR. Genes were selected on the basis of relationships with our drug targets, as shown in Fig. 4A, and included SST, CXCL8, NAMPT, TNF, NFKBIZ, BIRC3, TNFSF14, APOE, NKX3-1, HBA1, GPR182, SPDEF, STAB2, and DMBT1. Each of these genes either leads to increased activation of one or more of our drug targets or is transcriptionally regulated by one or more drug targets. Unexpectedly, the differential expression seen in our RNA-Seq was confirmed by qPCR validation in just three of the 14 selected genes (Fig. 4B). For two additional genes, there was significant differential expression between the primary tumor and liver metastasis, but not in the primary tumor and pooled metastasis comparison, while four genes displayed a trend toward higher expression in the liver metastases (0.05 < P < 0.1). Overall, nine of 14 genes trended toward or were significantly altered in pNET metastases versus primary tumors.

Compounds predicted by IPA and CMap to have suppressive activity against pNETs were tested using a quantitative, cell viability assay in two pNET cell lines, BON-1 and QGP-1. Fifteen drugs were selected for testing: triptolide, alvocidib (flavopiridol), mocetinostat, entinostat, bisindolylmaleimide-IX (RO 31–8220), PIK-75, PI-103, sirolimus (rapamycin), simvastatin, doxorubicin, daunorubicin, PD-184352 (CI-1040), PD-325901, selumetinib (AZD-6244), and apitolisib (GDC-0980). For each drug tested, the PNET cell lines displayed similar sensitivities in proliferation/viability assays although in most cases BON-1 cells were slightly more sensitive, as indicated by the lower IC₅₀s (Table 1; Fig. 5). All drugs significantly inhibited the proliferation and survival of PNET cells at increased concentrations although differential responsiveness to the various compounds was observed. Drugs with the highest activities had IC₅₀ values ranging between 1 and 80 nmol/L and included three different MEK inhibitors (PD-325901, selumetinib, and PD-184352), two TOP2A inhibitors (daunorubicin and doxorubicin), as well as triptolide, PIK-75, and alvocidib. In comparison, drugs with the lowest antiproliferative activity in PNETs had IC₅₀ values ranging between 130 nmol/L and 8–9 μmol/L, which included HDAC inhibitors (mocetinostat and entinostat), PI3K and mTOR inhibitors (apitolisib and PI-103), and the least potent drug tested, simvastatin, an HMG-CoA inhibitor. Nine other drugs with low CMap scores were tested, of which five showed very low activity against these cell lines, and four showed moderate activity (Supplementary Table S6).

The incidence of metastatic pNETs is increasing, and although there are a wide variety of treatment options for patients with these tumors, response rates and survival benefits have been modest (6, 8, 9). An improved understanding of the transcriptomic changes underlying the metastatic process may help to identify novel therapies to improve treatment outcomes. We present a complementary bioinformatics approach to identify novel therapies for pNETs. A previous study demonstrated the utility of a similar approach using pooled gastroenteropancreatic NET data and identified the HDAC inhibitor entinostat as a potential treatment for metastatic NETs (31). In our study, 15 compounds identified through our analysis pipeline demonstrated significant antiproliferative activity against pNET cell lines in vitro.

A critical step in identifying the metastatic expression profile involved filtering out normal liver and lymph node genes. We frequently observed that metastases to the liver and lymph nodes overexpress genes associated with the site of metastasis. For example, several of the most highly expressed genes in the liver metastases are related to normal hepatic metabolism or coagulation. Several explanations for this observation exist: (i) hormonal influence from the tumor microenvironment induces expression of these genes; (ii) activation of certain hepatic pathways is important for metastatic cells to colonize and survive in the liver; or (iii) these genes may result from contamination of normal hepatocytes in our tumor specimens. The tissue-specific expression was significantly more pronounced in the liver metastases than in the lymph nodes. This may reflect a stronger transcriptional “reprogramming” effect in the liver or likelihood that while tumors from all sites have some degree of lymphocytic infiltrate, hepatocytes would only be predicted to be found in liver metastases. Histologic review of the tumors in this study revealed that neuroendocrine cells were the most prevalent cell type, with tumor cells making up 70%–95% of the samples, and samples with less than 70% cellularity were excluded from the analysis. Nonetheless, lymphocytic infiltrate was seen in almost all samples, and projections of hepatocytes were occasionally observed to extend within the gross tumor margin. Regardless of the etiology of this tissue-specific expression, we were interested in defining gene expression changes common to metastases regardless of the site, which required filtering out tissue-specific genes. The preponderance of genes removed by the filter were normal hepatic genes, and their removal significantly influenced subsequent analysis in IPA and CMap.

The filtered, common metastatic expression profile included 626 differentially expressed genes when the significance threshold was set at P < 0.05 using nonparametric analysis. A number of software tools are available to analyze RNA-Seq data including DESeq2, voom and others (6, 23), each of which employs different statistical methods and makes different assumptions about the data to determine differential gene expression. For our analysis, differential gene expression was determined using the Wilcoxon rank-sum test, which makes fewer assumptions about the distribution of data and is relatively insensitive to outliers as compared with DESeq or voom.

Because of the large number of statistical tests performed in any RNA-Seq experiment, correction for multiple testing is typically performed to reduce the incidence of false positive results, resulting in an adjusted P value, also known as the corrected P value, FDR, or q-value. When such a correction was performed for our data, the calculated FDR for nearly all genes in our metastatic PNET expression profile exceeded 0.9. These values may reflect both the broad similarity in gene expression between primary tumors and metastases as well as within group heterogeneity in gene expression. Analysis using PCA and t-SNE showed that for most samples, the primary driver of differential expression was the tissue of origin. However, for a minority, clustering was primarily determined by characteristics specific to individual patients. As illustrated in Fig. 3B, there were six liver metastases that clustered with their corresponding primary and lymph node samples, rather than with the larger group of liver metastases. This observation, coupled with the relatively indistinct borders between the clusters, highlights the heterogeneous nature of the disease. Nonetheless, the overall expression pattern proved remarkably capable of predicting drugs with activity against pNET cell lines. Moreover, there was significant overlap in the drug classes and specific agents predicted to be efficacious using two separate prediction tools, the IPA analysis match and CMap.

To evaluate the reproducibility of the metastatic gene expression profile, we evaluated several genes by qRT-PCR in a different group of pNET samples. These genes were selected due to their relationships to several important cancer pathways and drug targets, including NF-kB, PI3K, mTOR, CDK, MEK, PKC, and HDAC1, as indicated by IPA. Confirmation of these specific genes, which were used to predict the activation states of drug targets, would help to bolster the predictions made based on their differential expression. Although most genes selected could be validated using qPCR, it should be noted that the fold changes in the RNA-Seq data were near twofold, which corresponds to approximately a single PCR cycle. Ideally, we would have tested more tumors, but had our validation set was limited by the number of tumors with primaries, nodal, and/or liver metastasis samples available that were not tested by RNA-Seq. Thus, due to the high calculated FDR and the failure of many genes to validate using qPCR, one should be cautious to not give undue weight to any one gene in the common metastatic profile. Instead, the overall patterns seen in the data provide more relevant, biologically useful information.

Compounds that inhibit a variety of distinct signaling pathways were identified by IPA and CMap as potential therapeutic options for patients with metastatic pNETs. Some of the drug targets identified were the usual suspects implicated in pNET tumorigenesis, including PI3K and mTOR (9, 17–19). Inhibitors of both kinases were the most highly represented drug classes with connectivity scores less than −90 in CMap and highly negative scores in the IPA analysis match, and the combined PI3K and DNA protein kinase (DNA PK) inhibitor PIK-75 was among the more efficacious drugs tested. This drug has shown anticancer activity in preclinical studies (32, 33). It is poorly soluble and has not progressed to clinical trials, but improvements in drug delivery may allow for human trials in the future (34). Other inhibitors of mTOR and/or PI3K that were tested in this study included apitolisib, PI-103, and rapamycin (sirolimus). Rapamycin had overlapping IC₅₀s with PIK-75 in the PNET cells but never achieved maximal efficacy, while apitolisib and PI-103 were about eight- to 10-fold less potent than PIK-75. Even so, apitolisib and rapamycin have both been used clinically. In a phase I trial, apitolisib demonstrated modest antitumor activity in patients with advanced solid tumors (35). However, in a phase II trial it was less effective than everolimus for the treatment of renal cell carcinoma, with a significantly higher rate of treatment-limiting adverse effects (36). Interestingly, while there were 12 mTOR inhibitors with connectivity scores below −90, the widely used everolimus was not among the drugs suggested by either platform, with a connectivity score of −15.18 on the CMap.

Other targets, such as HDACs and CDKs, have recently generated significant interest as potential targets for NET therapy: (31, 37) the HDAC inhibitor, entinostat, and the CDK inhibitor, ribociclib, are both currently being investigated in phase II trials in NETs (trial numbers NCT03211988 and NCT02420691, respectively). One CDK inhibitor, alvocidib, was tested in our study and was among the drugs that showed the highest antiproliferative activity. Although systemic toxicity is a concern, it has been shown to have activity against a wide variety of hematologic and solid malignancies in phase I and II trials (38). We tested two HDAC inhibitors, entinostat and mocetinostat, and found they had moderate antiproliferative activity in the pNET cell lines. Entinostat is currently being investigated in clinical trials for NETs; however, mocetinostat has not yet been studied in NETs. To date, phase II trials using mocetinostat for other malignancies, either as monotherapy or as part of a combined regimen, have demonstrated modest antitumor activity accompanied by high rates of adverse events (39–42).

Inhibitors of the MAPK pathway, which includes MAPK kinase (MEK) and the proto-oncogene BRAF, were also well represented among the compounds suggested by CMap and IPA. Three MEK inhibitors, PD-0325901, PD-184352, and selumetinib were selected for testing. All effectively inhibited proliferation of both pNET cell lines, with PD-0325901 having the lowest IC₅₀ of any drug tested. MEK inhibitors have not been previously used for the treatment of NETs; however aberrant activation of the MAPK pathway has been implicated in NET tumorigenesis (16). MEK inhibitors have shown excellent activity against BRAF- or NRAS-mutated melanoma as well as Ras-driven plexiform neurofibromas (43), and studies investigating their use in combination with targeted or immune-modulating therapy in melanoma and other solid tumors are ongoing (44). PD-184352, also known as CI-1040, was the first MEK inhibitor to enter clinical trials; however, its activity in humans was limited by poor solubility and rapid clearance (45). The antitumor activity of PD-0325901 and selumetinib has been supported in phase I and II trials including patients with melanoma and other solid tumors (46–48).

Four HMG-CoA reductase inhibitors, also known as statins, were suggested by CMap, although none were highly ranked by IPA, and simvastatin was selected for testing. The antiproliferative and antimetastatic effects of simvastatin have previously been demonstrated in QGP-1 and BON-1 (49), and our testing confirmed its activity against both cell lines. However, high concentrations in the μmol/L range were needed to achieve an antiproliferative effect.

Among the better acting anti-NET drugs identified in our study, triptolide is a naturally occurring compound that exerts its antitumor effects through a variety of mechanisms, including inhibition of NF-kB signaling and induction of apoptosis (50, 51). Although preclinical studies with the drug have been promising, triptolide is poorly soluble in water, limiting its clinical application. The development of a water-soluble analogue, minnelide (52), has allowed for clinical testing, and a phase II trial in pancreatic adenocarcinoma (NCT03117920) is currently in progress. Several inhibitors of DNA topoisomerase were suggested by CMap, and daunorubicin and doxorubicin were selected for testing. Although cytotoxic chemotherapy is less commonly used for well-differentiated NETs, doxorubicin in combination with streptozocin is among the standard treatment options for well-differentiated PNETs (53). Daunorubicin is another prototypical topoisomerase inhibitors that is commonly used as a first-line treatment for acute myeloid leukemia (54).

The final drug selected for testing was bisindolylmaleimide-IX. This drug is part of the bisindolylmaleimide class of PKC inhibitors and is described as belonging to this class within CMap. Unlike other drug classes, which are highly represented in our analysis (e.g., mTOR inhibitors), only two PKC inhibitors received scores below −90, one of which is bisindolylmaleimide-IX, and the other (CGP-60474) is also a CDK inhibitor. In addition, there were three PKC inhibitors with scores greater than 90, indicating that treatment with these compounds was associated with expression changes mirroring (rather than opposing) those seen in metastatic PNETs. This fact combined with the observation that many PKC isozymes may act as tumor suppressors rather than tumor promoters (55), suggests that bisindolylmaleimide-IX's antitumor activity may occur through a different mediator. Although it has not been noted in CMap, bisindolylmaleimide-IX has also been shown to inhibit DNA-topoisomerase (56). Given the number of topoisomerase inhibitors suggested by CMap, it is possible that bisindolylmaleimide-IX's predicted activity against NETs reflects its inhibition of topoisomerase, rather than PKC.

Fourteen of 15 drugs selected for testing in BON-1 and QGP-1 showed significant activity against both cell lines (IC₅₀ ≤ 1 μmol/L). For those compounds that have been studied in humans (entinostat, mocetinostat, alvocidib, apitolisib, PD-0325901, selumetinib, daunorubicin, doxorubicin), significant inhibition of proliferation was seen at concentrations that are readily attained in phase I or II trials. In comparison, the mTOR inhibitor, sirolimus (rapamycin), which has been used in numerous human trials plateaued at 60%–80% efficacy even at the highest doses tested; simvastatin only showed sensitivity when doses of 10 μmol/L were reached. Nine additional compounds were tested that had low CMap scores (Supplementary Table S6), five of which showed little to no activity and four showed at least moderate activity, including vinblastine (Vinca alkaloid), thapsigargin (Calcium ATPase), MLN-2238 (proteasome inhibitor), and rigosertib (PI3K and PLK1 inhibitor).

Hofving and colleagues tested 1,224 compounds (at 1 μmol/L concentration) from a commercially available inhibitory library against BON1 and QGP1 (30). They found a general sensitivity to MEK inhibitors, drugs targeting HSP90 and Aurora kinases, but not to HDAC inhibitors. Comparing these results to our RNA-Seq strategy is of interest. We achieved greater than 50% inhibition for both cell lines at ≤1 μmol/L concentration with 14 of the 15 (93%) drugs selected. Of these, Hofving and colleagues confirmed sensitivity of both BON1 and QGP1 cell lines (IC₅₀ at 1 μmol/L) to 3 of the 14 (PD032591, PI-103, and PIK-75), to one but not the other cell line in 4 (selumetinib, doxorubicin, PD-184352, PI-103), whereas five showed limited or no activity in both cell lines (daunorubicin, sirolimus, mocetinostat, entinostat, bisindolylmaleimide-IX); 2 of the drugs we found were not in their screen (triptolide, apitolisib). For the 12 drugs tested in common, their screen found sensitivity for 25% of those we selected, partial sensitivity (one cell line only) for 33%, and missed 42% of these potentially useful drugs. These results indicate that although these screens can be complementary, screening of these 1,224 compounds would have resulted in many of the drugs we found to effectively kill both cell lines being missed, and would have also markedly increased the number of potential false positives that would need further evaluation. Although these results are promising, further studies will be required to confirm the clinical utility of these drugs. As described above, the expression changes between primary and metastatic PNETs are frequently small in magnitude, and this combined with tumor heterogeneity leads to significant difficulty in finding statistically significant differential gene expression. Sequencing of additional samples and the use of techniques that account for differences in tumor cellularity, such as tissue deconvolution or single-cell RNA-Seq, may help to overcome this issue.

While the drugs tested showed activity against BON-1 and QGP-1 pNET cells, both lines have extremely high proliferative indices compared with more typical well-differentiated pNETs (57). A strength of this study is that we identified new NET drugs using gene expression data from patients with well-differentiated tumors rather than cell line–based profiling, but confirmation of drug activity against less proliferative pNET models will ultimately be required. We and others in the NET field are beginning to generate well-differentiated, patient-derived NET spheroid cultures (58), providing an ideal system for future drug testing. Once optimal acting drugs are confirmed in well-differentiated NET cell populations, further preclinical testing in a murine model of metastatic NETs would enable measurement of drug tolerability, antiproliferative effects, and antimetastatic activity in vivo.

In conclusion, we present a pNET metastatic gene signature that was determined from a comparison of patient-matched metastases and primary tumors. Bioinformatic analyses of the gene expression changes using IPA and CMap identified drugs with predicted activity against metastatic pNETs, which was verified by in vitro cytotoxicity assays in two pNET cell lines. Our methodology uses patient-derived expression data to inform and focus subsequent drug testing, providing mechanistic rationale for compound selection and effectively prioritizing drugs for preclinical screening prior to progression to clinical trials.

T.A. Braun is an employee/paid consultant for ImmortaGene. C. Chandrasekharan is an employee/paid consultant for Lexcion. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A.T. Scott, P.J. Breheny, P.H. Ear, B. Darbro, T.A. Braun, D.E. Quelle, J.R. Howe

Development of methodology: A.T. Scott, P.H. Ear, B. Darbro, S. Umesalma, C.K. Maharjan, J.R. Howe

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A.T. Scott, P.H. Ear, C.A. Kaemmer, C.K. Maharjan, C. Chandrasekharan, J.S. Dillon, T.M. O'Dorisio, J.R. Howe

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A.T. Scott, M. Weitz, P.J. Breheny, P.H. Ear, B. Darbro, B.J. Brown, T.A. Braun, C.K. Maharjan, D.E. Quelle, J.R. Howe

Writing, review, and/or revision of the manuscript: A.T. Scott, M. Weitz, P.J. Breheny, P.H. Ear, B. Darbro, T.A. Braun, S. Umesalma, D.E. Quelle, A.M. Bellizzi, C. Chandrasekharan, J.S. Dillon, J.R. Howe

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A.T. Scott, B.J. Brown, T.A. Braun, G. Li, J.R. Howe

Study supervision: D.E. Quelle, J.R. Howe

The authors thank Dr. David Gordon for his guidance using CompuSyn software for drug response analyses. This work was supported by the T32 grant CA148062 (to A.T. Scott) and Neuroendocrine Tumor SPORE grant P50 CA174521 (to B. Darbro, T.A. Braun, D.E. Quelle, A.M. Bellizzi, J.S. Dillon, T.M. O'Dorisio, and J.R. Howe). RNA-seq data presented herein were obtained at the Genomics Division of the Iowa Institute of Human Genetics which is supported, in part, by the University of Iowa Carver College of Medicine (Iowa City, IA) and the Holden Comprehensive Cancer Center (Iowa City, IA; NCI of the NIH under Award Number P30CA086862).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.