Introduction

Ameloblastoma is a slow-growing, locally invasive, benign epithelial odontogenic neoplasm. It is thought to be arise from SOX2-expressing dental lamina epithelium1, remnants of the tooth-forming enamel organ2. The tumor exhibits epithelial cells resembling pre-ameloblasts on a basement membrane in loosely arranged cells resembling stellate reticulum while the stroma consists of loose connective tissue. Although odontogenic tumors are relatively rare, they constitute 3.8% of head and neck pathology, of which 40–50% are ameloblastoma3,4. Occasionally, ameloblastomas show malignant features or transform into malignancy5 and in rare cases metastasize6. Current treatment modalities range from conservative enucleation to radical excision and vary according to tumor subtype and location7,8. High recurrence rates (50–80%) have been observed in cases of conservative treatment9. Consequently, and despite recent advances in imaging-assisted surgical margin localization, post-operative histological confirmation is still required. This forces surgeons to either act conservatively, risking the need for a second surgery, or act aggressively thus increasing morbidity10 and the need for extensive reconstructive surgery.

There are few genomics and transcriptomics studies of ameloblastoma, with most investigations focusing on candidate-genes. Moreover, different comparison tissues were used in the handful of microarray studies; including gingival tissue11,12, whole tooth buds13, dentigerous cysts14 and a universal human reference RNA15. Furthermore, whole tumor samples used in these studies included large portions of stromal tissue. In spite of the heterogeneity in comparison tissue, there have been advances in understanding tumorigenesis of ameloblastoma.

A recent study examining the whole transcriptome of ameloblastoma suggested the existence of distinct molecular subtypes12. It is envisaged that better understanding of the molecular basis of ameloblastoma can aid the identification of diagnostic and prognostic markers and may lead to the development of novel, personalized treatment protocols16. To address this knowledge gap, we embarked on this study aiming to characterize the transcriptome of neoplastic ameloblastoma tissue and identify relevant molecular pathways and genes, using whole genome microarray.

Methods

Tumor collection and preparation

The study was conducted in accordance with approved human subject research guidelines and was approved by the local institutional review board and the ethics committee of the University of North Carolina-Chapel Hill. Between 2005 and 2008, 2 fresh frozen samples were obtained during surgical resection and 15 formalin-fixed paraffin-embedded (FFPE) samples were retrieved from the archives of the Department of Oral and Maxillofacial Pathology Laboratory, University of North Carolina (UNC) School of Dentistry. All samples were evaluated by a board-certified oral and maxillofacial pathologist and at least one other author and diagnoses were classified based on the 2005 WHO Histologic Classification of Odontogenic Tumors. Additional demographic data including gender, age, race, and tumor recurrence were recorded and examined for potential associations.

The fresh tumors including the bony resected margins were placed in RNAlater and Richard Allan Scientifics’ decalcifying solution (water, hydrochloric acid, EDTA, tetrasodium tartrate and potassium tartrate) at 4 °C for 1–4 weeks. They were then frozen at −80 °C before 7 μm sections were obtained under RNAse-free conditions15, stained lightly with hematoxylin & eosin and sent for immediate laser capture microdissection (LCM). The FFPE blocks were decalcified before embedding in paraffin. The samples were sectioned, stained and sent for laser capture.

Laser capture microdissection

The AutoPixTM automated LCM system (Arcturus Engineering, Santa Clara, CA, USA) was used to isolate tumor cells (basal epithelial cells adjacent to the basement membrane). Images of the tissue sections including the captured regions were obtained before and after LCM.

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 TargetAmpTM 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 experiments17. 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 control18. 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) project21 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 analysis22 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.

Figure 1: The micrographs show the laser capture of the epithelial portion of an ameloblastoma sample.
figure 1

(A) – light micrograph of follicular ameloblastoma at 4X showing tumor epithelial follicles that are single-cell thick with surrounding stroma. Arrows points to tumor follicles, (B) – light micrograph of follicular ameloblastoma at 20X showing the laser capture outline of epithelial cells, (C) – remnants of the stroma tissue after LCM and (D) – captured cells on capsure cap. Scale bar: 100 μm.

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.

Figure 2: Cluster analysis used to determine reference tissue.
figure 2

(A) – Heat map of the 3 different odontogenic tissues (OB – Odontoblast, PA – Pre-secretory ameloblast, SA – Secretory ameloblast) and the 2 distinct clusters of Ameloblastoma (AM) clustered using the 60 odontogenic epithelium-defining genes. (B) – Array tree showing grouping of the 2 clusters with odontoblast and pre-secretory ameloblast. (C) – Participant demographic and tumor phenotype.

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 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).

Figure 3: Multiclass and pathway analysis of the different tumor clusters.
figure 3

(A) – Heat map of the genes with a FDR < 1% that are differentially expressed in the 4 clusters (OB – odontoblast, PA – pre-Secretory Ameloblast, oAM – odontoblast-like ameloblastoma, pAM – pre-secretory Ameloblast-like ameloblastoma) from a SAM multiclass analysis. The cluster tree on the right shows the 2 distinct clusters of ameloblastoma. Groups of genes with similar expression (identified by colored bars at the bottom of the heat map) were used for pathway analysis for the different clusters of ameloblastoma which were shown in (BD). (B) – Canonical pathways that are differentially expressed for the Common tumor cluster in IPA. (C) – Canonical pathways that are differentially expressed for the Pre-secretory ameloblast cluster in IPA. (D) – Canonical pathways that are differentially expressed for the odontoblast cluster in IPA.

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).

Table 1 Canonical pathways differentially expressed in ameloblastoma clusters.

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.

Table 2 Upstream analysis using Ingenuity Pathway Analysis.

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 populations23 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 development12. 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 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 mediators13. 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 extensively24. 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 DeVilliers15. Some investigators found the activation of beta-catenin downstream of Wnt signaling25 while others found that expression of various Wnt pathways differ among ameloblastoma with the canonical Wnt pathway being the main transduction pathway26. 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 expression11,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 gene27 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 transformation28 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 tumorigenesis29,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 ameloblastomas31. 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 ERK32. 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 results33. 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 reconstruction34.

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 RNA35 and affect the accuracy of microarrays. However, recent studies supported the use of such samples for gene expression analysis36 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.

Additional Information

How to cite this article: Hu, S. et al. Ameloblastoma Phenotypes Reflected in Distinct Transcriptome Profiles. Sci. Rep. 6, 30867; doi: 10.1038/srep30867 (2016).