Ameloblastoma Phenotypes Reflected in Distinct Transcriptome Profiles

Ameloblastoma is a locally invasive benign neoplasm derived from odontogenic epithelium and presents with diverse phenotypes yet to be characterized molecularly. High recurrence rates of 50–80% with conservative treatment in some sub-types warrants radical surgical resections resulting in high morbidity. The objective of the study was to characterize the transcriptome of ameloblastoma and identify relevant genes and molecular pathways using normal odontogenic tissue (human “dentome”) for comparison. Laser capture microdissection was used to obtain neoplastic epithelial tissue from 17 tumors which were examined using the Agilent 44 k whole genome microarray. Ameloblastoma separated into 2 distinct molecular clusters that were associated with pre-secretory ameloblast and odontoblast. Within the pre-secretory cluster, 9/10 of samples were of the follicular type while 6/7 of the samples in the odontoblast cluster were of the plexiform type (p < 0.05). Common pathways altered in both clusters included cell-cycle regulation, inflammatory and MAPkinase pathways, specifically known cancer-driving genes such as TP53 and members of the MAPkinase pathways. The pre-secretory ameloblast cluster exhibited higher activation of inflammatory pathways while the odontoblast cluster showed greater disturbances in transcription regulators. Our results are suggestive of underlying inter-tumor molecular heterogeneity of ameloblastoma sub-types and have implications for the use of tailored treatment.

RNA extraction and microarray. Laser-captured cells from each tumor were pooled and total RNA was isolated with the PicoPure RNA Isolation kit (Arcturus Bioscience, Santa Clara, CA, USA). The Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA) was used to assess the yield and quality of total RNA. Amplification was completed on all samples using TargetAmp TM 2-Round Aminoallyl-aRNA Amplification Kit (Epicentre Biotechnologies, Madison, WI, USA).
RNA was then analyzed using the Agilent Whole Human Genome Microarray 4 × 44 K G4112F (Agilent Technologies, Palo Alto, CA) containing 44 thousand 60-mer oligonucleotides representing over 41 thousand human gene transcripts. For this step, 200 ng of RNA was converted into labeled cRNA with nucleotides coupled to fluorescent dye Cy3 using the Low RNA Input Linear Amplification Kit (Agilent Technologies, Palo Alto, CA) according to the manufacturer's protocol. The Human Universal Reference RNA from Stratagene (Santa Clara, CA, USA) was coupled with Cy5. Cy3-labeled cRNA (1.65 ng) from each sample and the Cy5-labeled universal reference was hybridized to the Agilent whole genome array 41 k formatted chips. Data were extracted using Feature Extraction version 9.5 (Agilent Technologies, Palo Alto, CA). Background subtraction and Loess normalization were performed using default setting of the Agilent extractor. The use of the universal RNA facilitated the use of the dentome as a comparison. It acts as a technical intra/inter normalizing control, decreasing variability by measuring signal output ratio of experimental to reference RNA rather than relying on absolute signal intensity two-color hybridization experiments 17 . The dentome consists of odontogenic tissue (microdissected samples of human odontoblasts, pre-secretory ameloblasts and secretory ameloblasts) expression data from previous work that employed the universal reference as a normalizing control 18 . The data set included 4 samples of each type of odontogenic tissue. It was collected from 12 anterior tooth buds (incisor and canine) from 4 different fetuses with each fetus contributing 3 tooth buds. Each type of odontogenic tissue was collected from a single tooth bud with each of the 3 buds providing a single type of odontogenic tissue or developmental stage. (Gene Expression Omnibus microarray database accession number GSE63289). The ameloblastoma expression data were submitted to the Gene Expression Omnibus microarray database (accession number GSE68531).

Microarray data analysis.
A multiclass analysis was conducted between the 3 types of odontogenic tissue and the 60 genes differentially expressed at a false discovery rate (FDR) < 20% were designated as the odontogenic tissue-defining genes. The most appropriate comparison tissue was decided to be the normal tissue with the most similar profile to ameloblastoma, such that identified differences would be tumor specific. Cluster analysis was conducted using Cluster 3.0 between the 3 normal and tumor samples and visualized using Java TreeView-1.1.6r2.
Differential gene expression between tumors and comparison tissue were examined using Significance analysis of microarrays (SAM) 4.0. Ingenuity pathway analysis was used to identify differentially expressed pathways. In addition, upstream analysis from the ingenuity pathway analysis software was conducted. Gene set enrichment analysis (GSEA) 19 was conducted using GSEA v2.1.0 from the Broad institute (Cambridge, MA, USA) and the "all curated gene sets v4.0" available via the Molecular Signatures Database (MSigDB) 20 . Additionally, the ameloblastoma transcriptome was compared with the 13 cancer molecular subtypes from The Cancer Genome Atlas (TCGA) project 21 to investigate possible correlation with other known cancer types.
Microarray gene expression validation using NanoString. A variety of approaches have been used to validate microarray data in the literature and NanoString was selected for the study. A random subset of 3 ameloblastoma and 2 control odontogenic tissue samples was used to validate the microarray gene expression data. NanoString nCounter (Seattle, WA, USA) high throughput gene expression analysis 22 was performed using the Human Cancer Reference codeset (http://www.nanostring.com/products/gene_expression_panels). Each reaction contained 50 ng of total sample RNA plus reporter and capture probes. Digital counts were extracted, normalized and analyzed using nSolver v2.5 software. Differential expression between ameloblastoma and pre-secretory ameloblast from nanoString was compared with that obtained from the microarray.

Results
LCM facilitated the isolation of basal epithelial (neoplastic) cells from the tumor samples ( Fig. 1) without contamination from surrounding stroma cells. RNA was extracted from the LCM samples, with the 260/280 ratio for the 17 samples between 1.7-2.5 and a yield of between 88ng to 928ng (Supplemental Table 1). Quality control analysis conducted on the microarray chips indicated expected values for positive and negative controls, as well as uniformly high detected genes in both the red and green channels. Overall, no outliers were detected in either the normal tissue or tumor arrays indicating consistency in hybridization between samples.
NanoString validation. Due to the limited RNA quantity available, we did not conduct qPCR validation; instead, the NanoString platform was used to validate microarray gene expression data. Fold changes obtained with the nCounter system were correlated with those obtained from the microarray for the same samples (Supplementary Table 2). The 2 sets of expression data showed a good Pearson correlation (r = 0.61) in the scatter-plot ( Supplementary Fig. 1).
Determination of comparison tissue. Unsupervised hierarchical cluster analysis was conducted for the tumor and normal tissue samples which showed the presence of 2 distinct clusters of ameloblastoma and a separate cluster of normal odontogenic tissue ( Supplementary Fig. 2). Supervised cluster analysis using the 60 odontogenic tissue defining genes showed that the 2 clusters of ameloblastoma associated most closely with pre-secretory ameloblast (PA) and odontoblast (OB) ( Fig. 2A,B). As such, these 2 clusters were designated as the pre-secretory ameloblast-like ameloblastoma (pAM) and odontoblast-like ameloblastoma (oAM), respectively. Out of 10 samples in the pAM cluster, 9 were of the follicular type while 6/7 of the samples in the oAM cluster were of the plexiform type. (Fig. 2C). A Chi-Square analysis showed that the molecular clusters were significantly associated with a histological subtype (p < 0.05). A single comparison tissue could not be designated as the 2 clusters associated most closely with different odontogenic tissue; instead, a multiclass approach was employed.

Multiclass analysis.
Multiclass analysis was conducted between the 2 tumor clusters and 2 associated normal tissues (pAM, oAM, PA, OB) and differentially expressed genes were carried forward in a cluster analysis (Fig. 3A). To characterize the transcriptome of ameloblastoma and the 2 molecular sub-clusters, the gene  Continued expression data were analyzed in 3 groups. The common tumor cluster describes differential gene expression common in both tumor clusters compared to normal tissue and comprises 2592 genes that were expressed at a higher and lower level in the 2 tumor clusters (pAM, oAM) compared to the 2 normal tissue clusters (OB, PA). The pAM cluster describes differential gene expression unique to that cluster and consists of 1287 genes expressed at higher and lower levels compared to the other 3 groups. The oAM cluster describes differential gene expression unique to oAM tumors and consists of 1516 genes expressed at a higher and lower levels compared to the other 3 groups. The genes with fold changes at FDR < 1% were used for pathway analysis (Supplementary Table 3).
Pathway and gene set enrichment analysis. Ingenuity pathway analysis was used to examine the activated and inhibited canonical pathways for each tumor cluster ( Table 1). The common tumor cluster had 21 activated (z-score > 1) and 5 inhibited (z-score < − 1) pathways (Fig. 3B) at p-value < 0.05. Genes associated with notable biological processes that were differentially expressed in all the ameloblastoma tumors included prevention of damage to cell cycle regulation, cancer pathways, inflammatory pathways and Map kinase related pathways.
In addition, GSEA conducted between the common tumor cluster and normal tissues showed that 1860 out of the 2381 genes sets in the "all curated gene sets v4.0" were up-regulated in the common tumor cluster. Nineteen upregulated gene sets were significantly enriched at the nominal p-value < 0.05. (Supplementary Table 4).
The pre-secretory ameloblast tumor cluster had 22 activated and 1 inhibited pathway below the critical p-value threshold (Fig. 3C). Pathway analysis showed activation in the known cancer pathways, several inflammatory pathways and EGFR pathways.
The odontoblast tumor cluster had 1 activated and 8 inhibited pathways (Fig. 3D). Several inflammatory pathways were found to be inhibited in this cluster.
Upstream analysis. Upstream regulators that were predicted to be activated or inhibited are listed in Table 2. Most of the predicted upstream regulators in the common tumor cluster were transcription regulators, kinases and cytokines. Specifically, several Map kinase members and inflammatory cytokines were predicted to be activated.
Correlation with The Cancer Genome Atlas. The 2 molecular subtypes of ameloblastoma were compared with the transcriptome of the cancer subtypes in TCGA ( Supplementary Fig. 3). The analysis did not show any significant correlation of ameloblastoma with any of the 13 subtypes of cancers that are well studied and has established treatment protocols. As ameloblastoma does not seem to correlate molecularly with the cancer subtypes, more investigation into ameloblastoma tumorigenesis is needed for the development of effective treatment.

Discussion
The major finding of this study was the molecular heterogeneity of ameloblastoma that was strongly associated with its histological subtypes. Gene expression profiles of follicular and plexiform subtypes were more closely related to gene expression profiles of different normal odontogenic tissues and the follicular subtype showed activation of different molecular pathways compared with the plexiform subtype. This new knowledge can serve as a rich hypothesis-generating resource for the study of molecular and phenotypic characteristics of ameloblastoma.
One of the strengths of this study was the use of LCM for the examination of odontogenic tumors. The ability of LCM to isolate one cell-thick discrete tissue populations 23 facilitates the targeting and pooling of neoplastic epithelial portions of ameloblastoma. Similar to the present study, Heikinheimo and colleagues found that ameloblastoma gene expression is heterogeneous, and identified 2 distinct tumor clusters with gene expression profile that were most similar to gene expression in the cap/bell stage of tooth development 12 . Using supervised cluster analysis we found that more than half of ameloblastoma samples were most similar in gene expression to pre-secretory ameloblast, similar to those observed in the early cap/bell stage as described by Heikinheimo. The remaining ameloblastoma samples were associated with the mesenchymal-derived odontoblast rather than the epithelial-derived ameloblast; this appears to be driven by differences in inflammatory pathways and was associated with a different histological appearance. However, the early pre-secretory ameloblast in this sample could Canonical pathways different between the ameloblastoma common tumor cluster and normal cells  Table 1. Canonical pathways differentially expressed in ameloblastoma clusters. * p-value < 0.05.
be more appropriately described as preameloblast from the early cap stage of tooth development. Preameloblasts share a number of genes with odontoblasts due to significant cross-talk during differentiation; it is therefore unsurprising that a cluster of ameloblastoma was associated with odontoblast. The examination of pathways common to both tumors shows that inflammation appears to play an important role in ameloblastoma tumorigenesis and proliferation. This finding is similar to an earlier microarray study showing increased expression levels of inflammatory mediators 13 . Canonical pathway analysis showed that several immune/inflammatory pathways are activated in addition to several predicted activated upstream cytokines. The association between dysregulated inflammation and cancer progression has been studied extensively 24 . The greater number of pathways activated in the pre-secretory cluster suggests that this process may be more important in the follicular subtype.
The Wnt signaling pathway is important in the development of ameloblastoma as discussed in the microarray study by DeVilliers 15 . Some investigators found the activation of beta-catenin downstream of Wnt signaling 25 while others found that expression of various Wnt pathways differ among ameloblastoma with the canonical Wnt pathway being the main transduction pathway 26 . In this study, the canonical Wnt/β -catenin Signaling pathway was found to be activated in the common tumor cluster, highlighting its importance in tumorigenesis.
Additionally, both tumor clusters revealed that damage to cell cycle regulation pathways play important roles. Alteration in cell cycle regulation has been found by other investigators examining ameloblastoma gene expression 11,14 . A key regulator in the cell cycle damage prevention pathways is TP53 which is also predicted to be inhibited in our upstream analysis. TP53 is a major tumor suppression gene 27 and the loss of a tumor suppressor gene activity in ameloblastoma may be important in the tumorigenesis process. Although the dental epithelium defining SOX2 was not found to be differentially expressed, SOX11 of the same family of transcription factors related to oncogenic transformation 28 was found to be inhibited in the tumor samples.
Several canonical pathways involving the MAPK pathways and upstream members were found to be activated in the common tumor cluster. The MAPK pathways have long been considered tumor driver pathways in the pathogenesis of various cancers and are thought to be important in ameloblastoma tumorigenesis 29,30 . Although the expression levels of usual targets such as BRAF and EGFR were not directly increased, activation of the pathway suggests alterations in the activity of the members. Recently, Kurppa and colleagues found BRAF gene mutations, specifically V600E mutation, in 63% of ameloblastomas 31 . BRAF is in the RAS pathway and MAPK cascade. This mutation leads to a 500 fold increase in activity of BRAF and increases signals through MEK to activate ERK 32 . ERK is the activator of numerous downstream transcription factors which induces biochemical functions such as cell differentiation, proliferation, growth, and apoptosis. The V600E mutation in BRAF is a promising oncogene target for the anti-neoplastic drug dabrafenib; which was used in conjunction with a MEK inhibitor in a patient with stage 4 ameloblastoma with good results 33 . The use of chemotherapeutic agents to reduce tumor size can be very helpful in cases of ameloblastoma requiring extensive surgical margins and major post-surgical reconstruction 34 .
In GSEA, the gene sets with the highest activation scores, were cancer related including "SMID_ BREAST_CANCER_LUMINAL_A_DN". Cancer related pathways were also found to be differentially expressed in our pathway analysis. Moreover, breast cancer specific "Role of BRCA1 in DNA Damage Response" pathway was found to be differentially expressed, with the 2 different analyses highlighting similar molecular pathways.
One of the short-comings of this study is that most of the samples were FFPE. Formalin fixing can cause the degradation of RNA 35 and affect the accuracy of microarrays. However, recent studies supported the use of such samples for gene expression analysis 36 and NanoString has been shown to produce consistent results independent of the sample type (fresh frozen versus FFPE) 37 . Genes in our microarray data that had the greatest fold changes showed good correlation with the nanoString expression. In addition, there were no outliers in the cluster analysis among the ameloblastoma cluster analysis indicating consistent results between fresh frozen and FFPE samples.
Strengths of the study included the use of LCM and universal RNA as a means of normalizing between arrays. Ameloblastoma presents with neoplastic epithelial tissue surrounded by stromal tissue making isolation very difficult. As a result, most investigations of ameloblastoma used samples that contain diverse cell populations such as the surrounding stroma which can obscure driver pathways from the actual neoplastic epithelial cells. The use of the universal RNA facilitated the use of the dentome as a comparison and also allows the use of the microarray data by other investigators using universal RNA for normalization.
In conclusion, our study isolated ameloblastoma epithelial and normal odontogenic cells using LCM to identify gene expression profiles and molecular pathways that are potentially important in the tumorigenesis of ameloblastoma. Ameloblastoma showed 2 distinct molecular profiles that were associated with different histological  subtypes suggesting they could be receptive to different chemotherapeutic protocols. These results provide a wealth of information that can be used in future experimental and mechanistic studies, involving animal models and new pharmacogenomic approaches.