Characterizing and classifying neuroendocrine neoplasms through microRNA sequencing and data mining.

Neuroendocrine neoplasms (NENs) are clinically diverse and incompletely characterized cancers that are challenging to classify. MicroRNAs (miRNAs) are small regulatory RNAs that can be used to classify cancers. Recently, a morphology-based classification framework for evaluating NENs from different anatomical sites was proposed by experts, with the requirement of improved molecular data integration. Here, we compiled 378 miRNA expression profiles to examine NEN classification through comprehensive miRNA profiling and data mining. Following data preprocessing, our final study cohort included 221 NEN and 114 non-NEN samples, representing 15 NEN pathological types and 5 site-matched non-NEN control groups. Unsupervised hierarchical clustering of miRNA expression profiles clearly separated NENs from non-NENs. Comparative analyses showed that miR-375 and miR-7 expression is substantially higher in NEN cases than non-NEN controls. Correlation analyses showed that NENs from diverse anatomical sites have convergent miRNA expression programs, likely reflecting morphological and functional similarities. Using machine learning approaches, we identified 17 miRNAs to discriminate 15 NEN pathological types and subsequently constructed a multilayer classifier, correctly identifying 217 (98%) of 221 samples and overturning one histological diagnosis. Through our research, we have identified common and type-specific miRNA tissue markers and constructed an accurate miRNA-based classifier, advancing our understanding of NEN diversity.

INTRODUCTION

Classifying neuroendocrine neoplasms (NENs) is challenging due to tumor diversity, inconsistent terminology and piecemeal molecular characterization. Currently, NENs are broadly divided into epithelial or non-epithelial groups based on site of origin and differences in keratin and other gene expression; each group comprises multiple pathological types (1–3). To facilitate comparisons between NENs from different anatomical sites, international experts recently proposed a common classification framework (3). Here, the terms ‘category’, ‘family’, ‘type’ and ‘grade’, respectively, denote predominant neuroendocrine differentiation, degree of differentiation, diagnostic entity and inherent biological activity. While morphological assessment and immunohistochemical staining for chromogranin A, synaptophysin and Ki-67 proteins remain indispensable for confirming neuroendocrine differentiation and assessing tumor grade, other relevant molecular findings will be integrated into this framework over time. These studies will unravel many puzzles in NEN biology, including delineating the molecular differences between well-differentiated neuroendocrine tumors (NETs) and poorly differentiated neuroendocrine carcinomas (NECs) and finding regulatory molecules that underpin the ‘common neuroendocrine multigene program’ (3).

MicroRNAs (miRNAs) are small (19–24 nt) regulatory RNA molecules that can also be used to classify cancer (4,5). miRNAs are highly informative tissue markers because of their abundance, cell-type and disease-stage specificity, and stability in fresh and archived materials (6,7). These molecules also provide valuable mechanistic insights into cellular processes due to computationally predictable interactions with messenger RNAs (mRNAs) (8,9). In addition, miRNA expression profiles can be used to assess data reliability and to prioritize mRNA targets through further organization into miRNA cluster and sequence family datasets (10). To date, multiple miRNA profiling studies have been performed on single or limited combinations of NEN pathological types using different RNA isolation, detection and analysis methods (11). Although these differences complicate interstudy comparisons, miRNAs still hold much promise as multi-analyte markers that better reflect the ‘complexity and multidimensionality of the neoplastic process’ than current mono-analyte markers (12,13). Given recent advances in miRNA detection and analysis (14), we expect that substantial biological and clinically relevant insights into NEN biology will be gained through comprehensive miRNA profiling of multiple pathological types.

Through small RNA sequencing and data mining, we have generated reference miRNA expression profiles for multiple NEN pathological types and site-matched non-NEN controls, identified candidate category- and type-specific miRNAs, found evidence for constitutive and convergent miRNA gene expression in epithelial and non-epithelial NENs, and established a novel multilayer classifier for discriminating NEN pathological types.

MATERIALS AND METHODS

Study design and clinical materials

Sequencing-based miRNA expression profiles from 378 clinical samples, comprising 239 NEN cases and 139 site-matched non-NEN controls, were used in this study. Expression profiles were either compiled from published studies (7,15–18) (n = 149) or generated through small RNA sequencing (n = 229). Diagnostic histopathology, small RNA cDNA library preparation and the source of each sample are presented in Supplementary Table S1. The use of de-identified clinical data and banked or archived clinical materials was approved through the Research Ethics Board at Queen’s University, the Institutional Review Boards of Memorial Sloan Kettering Cancer Center, The Rockefeller University and Weill Cornell Medicine, and the Medical Ethics Committee at the Amsterdam University Medical Center.

RNA isolation and quantitation

Total RNA was isolated from 306 formalin-fixed paraffin-embedded tissue blocks and 72 fresh-frozen tissue samples using the Qiagen RNeasy® Mini Kit (n = 258), TRIzol™ Reagent (n = 68), the Ambion RecoverAll™ Total Nucleic Acid Isolation Kit (n = 28), Amsbio RNA-Bee™ Isolation Reagent (n = 10) and Qiagen miRNeasy® Mini Kit (n = 5), according to the manufacturers’ instructions or as described (7,15–18). Total RNA concentrations were measured using the Qubit™ fluorometer (n = 258), NanoDrop® ND-1000 spectrophotometer (n = 61) or Agilent 2100 Bioanalyzer (n = 28). RNA isolation and quantitation data were unavailable for 9 (2.4%) and 31 (8.2%) samples, respectively.

Small RNA sequencing and sequence annotation

miRNA expression profiles for all 378 samples were generated using an established small RNA sequencing approach and sequence annotation pipeline (10); spiked-in oligoribonucleotide calibrator markers enabled miRNA quantitation in each sample. Small RNA cDNA libraries were sequenced on HiSeq 2500 Illumina platforms in the Genomics Resource Center, The Rockefeller University, the McGill University and Génome Québec Innovation Center, and the Genomics Core, Albert Einstein College of Medicine. FASTQ sequence files were annotated through an automated pipeline (rnaworld.rockefeller.edu) (10), yielding sequencing statistics and merged miRNA, miRNA cluster and calibrator read counts. Merged miRNA refers to combined counts of multicopy miRNAs from different genomic locations and miRNA clusters are transcriptional units as defined (19); the term ‘miRNA’ will hereafter refer to merged miRNA data. Annotated sequencing statistics for all samples are presented in Supplementary Table S2; miRNA content was calculated using total RNA and calibrator RNA input ratio multiplied by total miRNA and calibrator count ratio as described (7). miRNA, miRNA cluster and calibrator read counts for all samples are presented in Supplementary Tables S3–S5, respectively.

Data preprocessing and filtering

Data preprocessing, filtering and subsequent analyses were performed in MATLAB (Mathworks, Inc., Natick, MA, USA, version R2019a) as described (18). To report miRNA abundance independent of sequencing depth, read counts were normalized against total sequence reads annotated as miRNAs. Sample outliers and batch effects were identified through correlation analyses (20) of miRNA expression profiles and excluded from the final dataset to increase study rigor. These analyses were completed for each NEN pathological type or site-matched non-NEN control group prior to preprocessing of the combined sample set. Sequencing data were of sufficient quality for 221 (92%) of 239 NEN cases and 114 (82%) of 139 non-NEN controls. Most excluded samples were individual outliers, except for 10 non-NEN samples from a single sequencing run. Following preprocessing, all non-human miRNAs and human miRNA STAR sequences were excluded from further analyses. To exclude miRNAs or miRNA clusters with low expression across samples, a filtering threshold was applied as described (6); specific filtering thresholds were set as a percentile of overall expression as indicated below.

Unsupervised hierarchical clustering of filtered miRNA expression profiles

To assess sample grouping, unsupervised hierarchical clustering was performed using log2 transformed normalized read counts of miRNA and miRNA clusters from all preprocessed samples. Euclidean distance was used as the similarity parameter with complete agglomeration clustering applied in the heatmap.2 function of the R gplots package (www.rdocumentation.org/packages/gplots/versions/3.0.1.1). Lowly expressed miRNAs and miRNA clusters were excluded with the filtering threshold set as the top 75% abundant miRNA and clusters in at least one sample.

Assessment and comparative analyses of abundant miRNAs in NEN and non-NEN samples

To identify candidate miRNA markers for all NENs and for each NEN pathological type, we ranked miRNAs and miRNA clusters by abundance and considered the top 1%. These abundant miRNAs and miRNA clusters were compared and correlated between NEN cases and non-NEN controls, as well as between each pathological type and site-matched non-NEN control group. To highlight substantial differences in miRNA expression, only comparisons with 20-fold or greater difference are discussed. For single-member miRNA clusters, abundance measures approximate the abundance of the single miRNA, and are not separately discussed.

Discovery analyses for miRNA-based NEN classification

To identify miRNAs or miRNA clusters that accurately discriminate between or within epithelial or non-epithelial NENs, we used an established feature selection algorithm that is an ensemble of 12 different machine learning techniques with 5-fold cross-validation (20). To prioritize high expression, we set the filtering threshold to the 90th percentile; miRNAs or miRNA clusters expressed above this threshold in >5% of samples were retained. Next, we ranked miRNAs and miRNA clusters that discriminate epithelial from non-epithelial NENs (comparison A). We subsequently ranked miRNA markers that successively identified epithelial NENs, including parathyroid adenoma (PTA), pituitary adenoma (PitNET), Merkel cell carcinoma (MCC), medullary thyroid carcinoma (MTC) and lung NENs from gastrointestinal–pancreatic (GEP) NENs (comparisons B–F), respectively. Lastly, we identified miRNA markers that discriminated neuroblastoma (NB), pheochromocytoma (PCC) and extra-adrenal paraganglioma (PGL) from each other (comparisons G and H) within the non-epithelial group. Only the top-ranking 3% miRNAs and miRNA clusters in these comparisons were assessed for classification below.

Construction and cross-validation of multilayer classifier

Scaling our existing approach to miRNA-based NEN classification (18,20), we constructed and cross-validated a multilayer classifier for discriminating NEN pathological types based on selected miRNAs. For each decision layer, all available algorithms (n = 23) in the MATLAB Classification Learner App were evaluated using 5-fold cross-validation. In each case, the classification model with highest accuracy was a support vector machine (SVM) classifier that was used to identify the smallest subset of miRNAs with the most discriminatory power for comparisons A–H above. Based on these subsets, we constructed a multilayer classifier through which miRNA profiles were first assigned as epithelial or non-epithelial prior to assignment to a specific pathological type.

Assessment of classifier performance and transferability

To assess the performance and transferability of our multilayer classifier, we used t-stochastic neighbor embedding (t-SNE) to visualize sample grouping patterns based on miRNAs selected for classification. We also determined overall classifier accuracy, evaluated the impact of miRNA cluster member substitutions on classifier accuracy and inspected the expression levels of the selected miRNAs.

Statistical analyses

Statistical analyses of clinical parameters were performed using SPSS Statistics (IBM, Armonk, NY, USA, version 25) and MATLAB. Differences in miRNA content and normalized miRNA expression were evaluated between NEN and non-NEN samples, and within NEN pathological types using the non-parametric Kruskal–Wallis (K–W) test (21); a two-sided P-value of <0.05 was considered significant. Similarities in miRNA expression between samples were determined using Spearman’s correlation (22).

RESULTS

Anatomical distribution and histopathological diagnoses of study samples

To characterize and compare miRNA expression between NEN and non-NEN samples, we collected relevant study materials, generated comprehensive miRNA expression profiles through barcoded small RNA sequencing, quality controlled profiles through data preprocessing and performed downstream analyses using statistical and machine learning approaches. Following data preprocessing for quality control, our final study cohort comprised 221 NEN cases and 114 site-matched non-NEN controls, hereafter termed study samples (Table 1). NEN cases comprised 15 distinct pathological types, arising in seven anatomical sites, including the gastrointestinal tract and pancreas, lung, parathyroid gland, pituitary gland, skin, thyroid gland, and the adrenal gland and extra-adrenal sites. Site-matched non-NEN controls comprised non-diseased tissues and non-NEN cancers from five anatomical sites, including the pancreas, lung, parathyroid gland, skin and thyroid gland.

Table 1.

Anatomical distribution and histopathological diagnoses of study samples

NENs	Number of samples, n (%)	non-NENs	Number of samples, n (%)
Total	221		114
Epithelial
Gastrointestinal
tract and pancreas
PanNET	28 (13%)	PAAD	10 (9%)
INET	31 (14%)
AppNET	15 (7%)
RNET	7 (3%)
Lung
TC	13 (6%)	LAC	9 (8%)a
AC	15 (7%)	LUNG	15 (12%)a
SCLC	11 (5%)
LCNEC	13 (6%)
Parathyroid gland
PTA	9 (4%)	PTG	15 (13%)
Pituitary gland
PitNET	10 (5%)
Skin
MCC	17 (8%)	SK	10 (9%)
Thyroid
MTC	9 (4%)	TG	10 (9%)b
		TN	45 (39%)b
Non-epithelial
Adrenal gland and
extra-adrenal sites
NB	25 (11%)
PCC	10 (5%)
PGL	8 (4%)

aFor lung NENs, neoplastic (LAC) and non-diseased (LUNG) tissue controls were available.

bFor MTC, neoplastic (TN) and non-diseased (TG) tissue controls were available.

Anatomical location and diagnostic histopathological information are presented for 221 NEN cases, comprising 15 pathological types from seven anatomical sites, and 114 site-matched non-NEN controls, comprising seven diagnostic entities from five anatomical sites. Sample abbreviations: AC, atypical carcinoid; AppNET, appendiceal NET; INET, ileal NET; LCNEC, large-cell NEC; MCC, Merkel cell carcinoma; MTC, medullary thyroid carcinoma; NB, neuroblastoma; PanNET, pancreatic NET; PCC, pheochromocytoma; PGL, paraganglioma; PitNET, pituitary adenoma; PTA, parathyroid adenoma; RNET, rectal NET; SCLC, small-cell lung carcinoma; TC, typical carcinoid. Non-NEN samples comprise lung (LUNG), lung adenocarcinoma (LAC), pancreatic adenocarcinoma (PAAD), parathyroid gland (PTG), skin (SK), thyroid gland (TG) and thyroid neoplasm (TN).

Small RNA sequencing of study samples

We generated comprehensive miRNA expression profiles for all samples through barcoded small RNA sequencing. Following sequence annotation, we obtained a median of 4 386 727 (range: 53 516–40 305 4453) total small RNA reads and 258 932 (range: 1312–3 723 507) calibrator reads (Supplementary Table S2). For miRNAs, we detected a median of 2 322 722 (range: 1740–34 781 174) miRNA sequence reads, representing a median of 63.1% total sequence reads; miRNA, miRNA cluster and calibrator expression profiles for each sample were subsequently generated from these reads. Median miRNA content was 26.4 (range: 0.0–2048.4) fmol/μg total RNA (Supplementary Table S2).

Abundant miRNA composition in NEN and non-NEN samples

To better understand miRNA composition in NEN and non-NEN samples, we assessed and correlated abundant miRNAs and miRNA clusters within and between sample sets. Abundant miRNA and miRNA cluster composition was similar within all NEN cases or all non-NEN controls. The number of members in each miRNA cluster is indicated in parentheses following the cluster name, e.g. cluster-hsa-mir-98(13). Among all NEN cases, miR-375, -21, -143, -let-7a, -26a, -7, -let-7f and -125b and cluster-mir-375(1), -98(13), -21(1) and -23a(6) were the most abundant miRNAs and miRNA clusters, with the median relative frequency ranging 1.5–10.6% and 3.6–10.6% of respective total read counts (Supplementary Table S6). Within this group, miR-375, -21, -143, -let-7a, -26a, -7, -let-7f, -125b and -141 and cluster-mir-98(13), -mir-375(1), -mir-7-1(3) and -mir-143(2) were highly expressed in five or more pathological types (Supplementary Table S6). In comparison, among all non-NEN controls, miR-21, -125b, -let-7a, -143, -let-7f, -30a, -26a and -29a and cluster-mir-98(13), -21(1), -30a(4) and -23a(6) were the most abundant miRNA and miRNA clusters, ranging 2.5–10% and 5.2–15.9% of respective total miRNA-annotated read counts (Supplementary Table S7). Within this group, miR-21, -let-7a, -143, -30a, -let-7b and -30d and cluster-mir-98(13), -mir-21(1), -mir-23a(6) and -mir-30a(4) were highly expressed in five or more non-NEN entities (Supplementary Table S7). Correlation analyses highlighted the similarities in abundant miRNA composition within epithelial and non-epithelial NENs; with the exception of PTA, NEN cases were less correlated with site-matched non-NEN controls (Supplementary Figure S1).

Comparative analyses of abundant miRNAs in NEN and non-NEN samples

To better understand meaningful differences in miRNA composition between NEN and non-NEN samples, we compared abundant miRNAs and miRNA clusters for all NEN samples and for each pathological type with relevant controls. Comparative analyses indicated that miR-375 and miR-7 were 216- and 48-fold higher in all NEN cases compared to all non-NEN controls, respectively. Fold changes ranged 59–816- and 41–69-fold higher for miR-375 and miR-7 in specific NEN pathological types compared to site-matched non-NEN controls (Supplementary Table S6 and Figure 1). The only exception was observed in PTA, which showed the lowest miR-375 and miR-7 expression of all NENs; in fact, higher expression was observed in non-neoplastic parathyroid glands. Other notable miRNA overexpression among NENs included miR-127, with 86-fold higher expression in typical carcinoids (TC) compared to lung non-NEN tissues (Supplementary Table S6); cluster-mir-127(8) was also 78-fold higher in TC compared to lung non-NEN tissues (data not shown). In addition, miR-203 and miR-205 expression was 143- and 366-fold higher in non-NEN skin controls than MCC, respectively (Supplementary Table S7).

Figure 1.

miR-375 and miR-7 expression in NEN and non-NEN samples. Normalized miR-375 and miR-7 expression was examined between 15 NEN pathological types and 7 site-matched non-NEN control groups. Site-matched NEN and non-NEN groups were designated by anatomical site in the color bar: pancreas (blue), lung (red), parathyroid (purple), skin (orange) and thyroid (green); NENs without a site-matched control were left blank. Both miR-375 and miR-7 were higher expressed in NEN cases than non-NEN controls. With the exception of PTA, miR-375 expression was higher in NEN pathological types than in site-matched non-NEN controls. With the exception of PTA, miR-7 was also higher in NEN pathological types compared to site-matched non-NEN controls. Abbreviation: log2 RF, log2 normalized relative frequency. Sample abbreviations are provided in Table 1 and Supplementary Table S1.

To assess the classificatory potential of miRNA expression profiling, we first explored our data using unsupervised hierarchical clustering. With the exception of all PTA samples and one large-cell NEC (LCNEC) sample, NEN cases and non-NEN controls clustered separately (Figure 2). In addition, epithelial samples clustered distinctly from non-epithelial samples with the exception of one pancreatic NET (PanNET). NEN pathological types preferentially clustered together rather than with site-matched non-NEN controls. Unsupervised hierarchical clustering of filtered miRNA cluster expression from the same samples clustered as above (Supplementary Figure S2).

Figure 2.

Unsupervised hierarchical clustering of study samples based on miRNA expression. Unsupervised hierarchical clustering using Euclidean distance and complete agglomeration clustering was performed using filtered (union of top 75% abundance) log2 normalized miRNA sequence reads for all NEN cases (n = 221) and non-NEN controls (n = 114). Anatomical groupings comprise the following pathological types described in Table 1 and Supplementary Table S1: thyroid (MTC, TG, TN), skin (MCC, SK), pituitary gland (PitNET), parathyroid gland (PTA, PTG), lung (AC, TC, SCLC, LCNEC, LAC, LUNG), GEP (AppNET, INET, PNET, RNET), and adrenal and extra-adrenal (PCC, PGL). With noted exceptions, NEN cases and non-NEN controls, and epithelial and non-epithelial samples, clustered distinctly and NEN pathological types preferentially clustered with each other than with site-matched non-NEN controls.

Next, we identified candidate miRNA markers for NEN classification using an established approach comprising feature selection and validation (18). Using this approach, we selected effective miRNA markers from the top-ranked 3% miRNAs or miRNA clusters discriminating between or within epithelial or non-epithelial NENs (Supplementary Tables S8 and S9). These comparisons were used to construct and assess the reliability of the multilayer classifier below.

We subsequently constructed and assessed the accuracy of a multilayer miRNA-based classifier for predicting NEN pathological types with 5-fold cross-validation (Figure 3). The resulting classifier consisted of eight decision layers, using the linear or cubic SVM model at each layer (Supplementary Table S10). In the first layer, miR-200a expression was significantly higher in epithelial than non-epithelial NENs (K–W P-value <0.01); miR-10b provided additional prediction power (K–W P-value <0.01). When combined, these two miRNAs discriminated epithelial and non-epithelial NENs with only one sample misclassification (Figure 4A), which was found to be a histological misidentification (see below). In subsequent layers, sample profiles were successively assigned to other pathological types using the least number of miRNAs required. Within epithelial NENs, PTA, PitNET, MCC and MTC were, respectively, discriminated from remaining NENs based on expression of miR-30a, miR-10a and miR-212-3p, miR-15b and miR-660, and miR-335-5p, miR-29a and miR-222 (Figure 4B–E). Lung NENs and GEP NENs were discriminated based on expression of miR-760, miR-1224-5p, miR-139, miR-205 and miR-9 with three misclassifications (Figure 4F and G). Within non-epithelial NENs, NB and PCC/PGL, and PCC or PGL, were accurately discriminated based on expression of miR-93, and miR-10b and miR-397, respectively (Figure 4H and I). Decision node level accuracy ranged from 97% to 100% (Supplementary Table S10).

Figure 3.

Multilayer miRNA-based classifier for predicting NEN pathological types. A multilayer classifier for predicting NEN pathological types was developed using supervised machine learning models. In the first layer, NEN miRNA profiles were classified as epithelial or non-epithelial based on miR-10b and miR-200a expression. In subsequent layers, epithelial and non-epithelial NENs were successively identified using the selected miRNAs as indicated. Sample abbreviations are provided in Table 1 and Supplementary Table S1.

Figure 4.

Scatter plot assessment of miRNAs selected for classification. Epithelial and non-epithelial NENs are effectively discriminated based on miR-10b and miR-200a expression with one misclassification (A). Within epithelial NENs, PTA, PitNET, MCC and MTC were accurately discriminated from the remaining NENs based on miR-30a expression (B), miR-10a and miR-212-3p expression (C), miR-15b and miR-660 expression (D), and miR-335-5p, miR-29a and miR-222 expression (E); lung NENs and GEP NENs were discriminated based on miR-760, miR-1224-5p, miR-139, miR-205 and miR-9 expression (F, G). Within non-epithelial NENs, NB was accurately discriminated from PCC/PGL based on miR-93 expression (H), and PCC and PGL were separated based on miR-10b and miR-379 expression (I). Similar results were generated using relevant miRNA cluster data and are not presented. Arrows indicate misclassified samples. Abbreviation: log2 RF, log2 normalized relative frequency. Sample abbreviations are provided in Table 1 and Supplementary Table S1.

Using the 17 miRNAs selected for multilayer classification, t-SNE analysis indicated clear separation of epithelial and non-epithelial NENs with one notable exception (Figure 5), which was found to be a histological misidentification (see below). NEN pathological types also grouped together within epithelial and non-epithelial clusters. With 217 of 221 samples accurately classified, the overall accuracy of our multilayer classifier was 98% (Table 2). miRNA cluster substitutions had little to no effect on overall and decision node level accuracy (data not presented). At each decision node of the classifier, selected miRNAs were always more highly expressed in one comparison group (0.40%; range: 0.01–7.82%) versus the other (0.03%; range: 0.00–2.35%; Supplementary Table S11), highlighting their potential as translatable tissue markers of specific NEN pathological types.

Figure 5.

t-SNE for selected classificatory miRNAs. Sample grouping was visually assessed using miRNAs selected for multilayer classification and t-SNE analysis. With one notable exception, samples clustered as epithelial or non-epithelial NENs and tended to group by pathological type. The exception was a misdiagnosed PanNET later found to be a PGL on further testing. Sample abbreviations are provided in Table 1 and Supplementary Table S1.

Table 2.

Overall accuracy of multilayer classifier for discriminating NENs

		Established diagnosis
		GEP NET	Lung NET	MTC	MCC	PitNET	PTA	PCC	PGL	NB
Multilayer classifier designation	GEP NET	80	3
	Lung NET		49
	MTC			9
	MCC				17
	PitNET					10
	PTA						9
	PCC							10
	PGL	1							8
	NB									25
Decision-level accuracy		80/81 (99%)	49/52 (94%)	9/9 (100%)	17/17 (100%)	10/10 (100%)	9/9 (100%)	10/10 (100%)	8/8 (100%)	25/25 (100%)
Overall accuracy		217/221 (98%)

Using our multilayer classifier, NEN miRNA profiles were assigned to one of nine pathological subgroups or pathological types. Cases of GEP NENs (AppNET, INET, PNET, RNET) or lung NENs (TC, AC, SCLC, LCNEC) were not assigned to individual pathological types because we previously developed miRNA-based classifiers for these subgroups (18) (Wong et al., in preparation). By comparing classifier designations to established histopathological diagnoses, we determined our overall classifier accuracy to be 98%. Additional measures of classifier performance were also calculated: precision (0.98), recall (0.99) and Matthews correlation coefficient (0.98). Sample abbreviations are provided in Table 1 and Supplementary Table S1.

Detection of histological misidentification by miRNA-based NEN classifier

The unusual finding of an epithelial PanNET within the cluster of non-epithelial NENs (Figures 4A and 5), in addition to miRNA-based classification of this PanNET as a PGL (Table 2), prompted us to review the histopathology of this case. Upon review, the tumor was a small (<1 cm in size) low-grade NET at the tail of the pancreas, with histological features overlapping both PanNET and PGL. Immunohistochemical analysis showed that the tumor cells were diffusely positive for synaptophysin, chromogranin A and GATA3, and negative for cytokeratin (AE1/AE3 antibodies). This phenotype diagnosed this tumor as a PGL, as predicted by the miRNA classifier, and not a PanNET, which should be cytokeratin-positive and GATA3-negative (23,24). The unusual case was misidentified based on initial histology, but was correctly diagnosed by molecular profiling and miRNA-based classification.

DISCUSSION

Accurate NEN classification is essential for understanding tumor biology and guiding clinical care. NEN pathological classification is modified by experts on an ongoing basis as updated clinical, pathological, biological and molecular data become available. Recently, these experts proposed a common classification framework for evaluating NENs, clarifying terminology to reduce confusion and harmonizing concepts to facilitate comparisons between pathological types (3). Although morphology-based, this framework is designed to incorporate ‘equally solid genetic studies across all anatomical sites (3)’ over time. Here, we generate biological and clinical insights into NENs through miRNA-based classification.

The strength of our study stems from comparing multiple NEN pathological types and site-matched non-NEN controls using comprehensive miRNA detection through barcoded small RNA cDNA library sequencing (25) and accurate sequence annotation (19). Advanced computational approaches for miRNA feature selection (20) and classifier construction (18) further bolstered our approach. We carefully assessed data reliability through knowledge of miRNA cluster composition (10), evaluated classifier performance and transferability by determining overall and decision node level accuracy, gauged the impact of miRNA cluster substitutions on accuracy and inspected the abundance of selected classificatory miRNAs. Throughout the study, miRNA clusters measured data quality and transferability of miRNAs as clinical markers; we then focused on miRNAs to build a streamlined prototype of a tool for NEN classification. The identified miRNAs can be used as mono-analyte or multi-analyte markers as needed (12).

Unsupervised hierarchical clustering of filtered miRNA expression profiles confirms existing knowledge and provides new knowledge of NEN grouping. With the exception of all PTA and one LCNEC sample, NEN cases and non-NEN controls clustered separately. Based on these findings, we speculate that all PTA have a distinct gene expression pattern linked to their indolent behavior; the LCNEC sample showed areas of possible squamous cell differentiation (data not shown) that may explain peculiar clustering patterns. Within NENs, two major groups corresponding to epithelial and non-epithelial NENs are evident; interestingly, one epithelial NEN clusters with non-epithelial NEN samples. Here, we show that these epithelial and non-epithelial NENs can be discriminated through miR-200a (26) and miR-10b expression, and confirm that our epithelial PanNET sample is actually a non-epithelial PGL based on additional cytokeratin and GATA-3 immunostaining (23,24). Within non-NENs, samples group mostly by anatomical site of origin as expected (6). Visual inspection of cluster diagrams indicates similarities and differences in abundant miRNA composition in NEN and non-NEN samples.

Similarities in abundant miRNA composition between samples provide coarse insights into cellular gene expression programs. Within NENs, miR-375, -21, -143, -let-7a, -26a, -7, -let-7f, -125b and -141 were highly expressed in five or more pathological types; known oncogenic or tumor suppressor functions for these miRNAs are reviewed elsewhere (8,27). miR-375, the most abundant miRNA in NENs, is believed to regulate lineage-specific differentiation (28–31), growth (32,33) and function (32,34) of neuroendocrine cells. Correlation analyses highlighted similarities in abundant miRNA composition for all NENs, including epithelial or non-epithelial NENs. These findings indicate that all NENs have a constitutive miRNA gene expression program that likely directly or indirectly maintains the neuroendocrine cell phenotype. Given the different cellular origins of epithelial and non-epithelial NENs (35), convergent miRNA gene expression likely implies functional similarities. Within non-NEN samples, miR-21, -let-7a, -143, -30a, -let-7b and -30d were highly expressed in five or more non-NEN entities; their cancer-related functions are reviewed elsewhere as above. While mechanistic insights into cellular processes can be gained through predictable targeting of mRNAs by abundant miRNAs, this topic is beyond the scope of the present study (36).

Differences in abundant miRNA composition between samples can also be used to identify new and confirm known miRNA markers. miR-375 expression was substantially higher in all NEN cases compared to non-NEN controls. Where comparisons allowed, miR-375 was consistently higher in NEN pathological types compared to site-matched non-NEN controls. Based on current miRNA expression tissue atlases, miR-375 is currently thought to be an endocrine gland specific marker (6,37). However, the presence of miR-375 in enteroendocrine cells (30,38), pancreatic beta cells (32,33), thyroid C cells (39), and MCC (7,31,40), NB (15) and SCLC cell lines (29) suggests that miR-375 is a neuroendocrine cell marker. Given the specificity and distribution of miR-375 in our samples and its reported abundance in seemingly disparate NEN pathological types (7,18,38,41–43), we propose that miR-375 is a universal marker of neuroendocrine cell differentiation. miR-375 appears to be highly expressed in NENs, in amounts proportional to the number of normal neuroendocrine cells and/or the degree of neuroendocrine differentiation within control tissues or tumors; neuroendocrine differentiation of tumors is more common than currently appreciated (44). More systematic studies are required to confirm this proposal.

Although less abundant than miR-375, miR-7 expression was also elevated in all NEN cases compared to non-NEN controls. Where comparisons allowed, miR-7 was often higher in NEN pathological types compared to site-matched non-NEN controls. Other than expression in the pituitary gland, atlas studies provide limited information on miR-7 expression (6,37). However, the presence of miR-7 in enteroendocrine cells (30), pancreatic islet cells (33,45), thyroid C cells (46), but not controls suggests that this miRNA also has some degree of neuroendocrine specificity. Given their specificity, some tissue profiling studies may have inadvertently interpreted miR-375 or miR-7 reduction in expansile cancer lesions as miRNA reduction rather than neuroendocrine cell destruction. Although miR-127 was higher in TC than non-NEN controls, the significance of this difference is unclear. Conversely, comparisons of abundant miRNA composition between non-NEN and NEN samples identified known tissue-specific miRNA markers such as miR-203 and miR-205 (6).

As with other cancers (4), miRNAs can be used for NEN classification. Using our feature selection algorithm, we identified 17 miRNAs to discriminate 15 NEN pathological types; t-SNE analyses using these miRNAs clearly separated epithelial and non-epithelial NENs and suggested clustering by pathological type. Given their classificatory potential, we subsequently constructed and validated a multilayer classifier for discriminating NEN pathological types, correctly identifying 217 (98%) of 221 samples. Three of the four misclassified samples occurred at the GEP NEN versus lung NEN decision node, suggesting model overfitting and the need for additional samples for validation. On further testing, the fourth ‘misclassified’ sample turned out to be a PGL as indicated by miRNA expression profiling. We also introduced criteria for evaluating classifier performance and transferability, including determining overall and decision node level accuracy, assessing the impact of miRNA cluster substitutions on classifier accuracy and showing the relative abundance of miRNAs selected for classification.

This study does have limitations that are commonly encountered in rare cancer and miRNA research. Comprehensive clinical information is challenging to obtain, limited sample numbers preclude hold out validation and miRNA content measurements can vary widely due to technical challenges. Nonetheless, we provide compelling evidence that miRNAs are useful for NEN classification and should be included in further multi-omic studies of these neoplasms.

Through comprehensive miRNA expression profiling, we have identified candidate universal and classificatory markers that may be useful as adjunct tissue markers, constructed a multilayer classifier for discriminating NENs and provided reference profiles for hypothesis generation or interstudy comparisons. Our next steps involve confirming our findings in well-annotated sample sets, evaluating miRNAs as circulating markers and investigating upstream promoter activity and downstream targeting events.

DATA AVAILABILITY

Annotated miRNA and miRNA cluster counts (Supplementary Tables S3 and S4) have also been deposited to Data Dryad (https://datadryad.org/stash/dataset/doi:10.5061/dryad.fn2z34tqj).

Click here for additional data file.