Microarray and bioinformatic analysis of conventional ameloblastoma: an observational analysis

Abstract

Ameloblastoma is a highly aggressive odontogenic tumor, and its pathogenesis is associated with many participating genes.

Objective

We aimed to identify and validate new critical genes of conventional ameloblastoma using microarray and bioinformatics analysis.

Methodology

Gene expression microarray and bioinformatic analysis were performed using CHIP H10KA and DAVID software for enrichment. Protein-protein interactions (PPI) were visualized using STRING-Cytoscape with MCODE plugin, followed by Kaplan-Meier and GEPIA analyses that were used for the candidate’s postulation. RT-qPCR and IHC assays were performed to validate the bioinformatic approach.

Results

376 upregulated genes were identified. PPI analysis revealed 14 genes that were validated by Kaplan-Meier and GEPIA resulting in PDGFA and IL2RA as candidate genes. The RT-qPCR analysis confirmed their intense expression. Immunohistochemistry analysis showed that PDGFA expression is parenchyma located.

Conclusion

With bioinformatics methods, we can identify upregulated genes in conventional ameloblastoma, and with RT-qPCR and immunoexpression analysis validate that PDGFA could be a more specific and localized therapeutic target.

Ameloblastoma; Computational Biology; Platelet-derived growth factor Alpha; IL2RA protein, human

Introduction

Odontogenic tumors (OTs) are oral lesions that impact the quality of life of patients because they affect not only the teeth but the maxilla and mandible. OTs constitute a group of heterogeneous diseases ranging from hamartomatous lesions to benign and malignant neoplasms with metastatic potential. These are derived from the epithelium, ectomesenchyme, and/or mesenchymal elements of the odontogenesis apparatus.¹

The epidemiology of OTs varies throughout the world; in some countries, the most frequent OT is ameloblastoma (Hong Kong, Japan, Zimbabwe, and Nigeria), whereas in others (United States of America, Brazil, and Canada) the most frequent tumor is odontoma.¹,² The most common OT in Mexico is odontoma, followed by ameloblastoma, myxoma, adenomatoid odontogenic tumors, and calcifying odontogenic cysts.³

Ameloblastoma is a slow-growing locally invasive benign OT with different histological variants, which may be located in the posterior zone of the mandible or maxilla. Because of its potential for recurrence, it is often classified as an aggressive tumor. The estimated global incidence of ameloblastoma is 0.5 cases per million people per year, with the age of diagnosis ranging from 30–60 years. Conventional ameloblastoma (CAm) is the most common variant, followed by unicystic and peripheral ameloblastoma. Although the precise etiology of CAm is unknown, dysregulation of many genes associated with odontogenesis is speculated to play an important role in its histogenesis.²,⁴ Changes in the expression or mutations in genes, such as BRAF, Ras, FGFR2, and SMO, among others, could be associated with its histogenesis.⁵,⁶ Given this complexity, high-throughput assays offer an alternative to comprehensively analyze this neoplasm. Microarray technology has been used to obtain information on the genetic alterations that occur in several diseases, including neoplasms, such as CAm. Much data are obtained with high-throughput analysis and integrated bioinformatics methods are necessary to unravel the mechanisms underlying the pathogenesis of diseases and to explore and identify novel biomarkers that could help us in further studies.⁷,⁸ Previous approaches with microarrays assays had shown elements of SHH, cell-cycle regulation, inflammation, MAP kinase pathways, and other molecules, which were confirmed via tests, such as PCR, immunohistochemistry, or NanoString, suggesting that these regulators are important elements of the pathogenesis of conventional ameloblastoma.⁹-¹¹ The bioinformatic analysis is only the first step for new biomarkers to be proposed. Those must be corroborated with particular assays, so that this information can cross over to the clinical level.¹² The objective was to identify and validate new critical genes of conventional ameloblastoma using microarray and bioinformatics analysis.

Methodology

Selection of CAm cases

This study was approved by the Institutional Technical Committee of the Support Program for Research and Technological Innovation Projects (DGAP/1956/2019) of the National Autonomous University of Mexico (UNAM). CAm samples were retrieved from the histopathological paraffin block archives of the Oral Medicine and Pathology Department, Postgraduate Division of Dentistry School (ISO-9001:2015 certified CMX C SGC 157 2017). This study was conducted following the integral privacy notice for patients from the Dentistry School, protecting their identity.¹³ All formol fixed paraffin embedded samples derived from patients who provided an informed consent form. The diagnosis was confirmed based on the 2017 World Health Organization (WHO) histological classification of OTs.² A total of 15 CAm and 16 dental follicles (used as controls) samples were obtained (Supplementary Table 1).

RNA extraction and cDNA microarray

Tissues (50 μm) were obtained from each sample.¹⁴ RNA extraction was performed using the ReliaPrep™ FFPE Total RNA Miniprep System (Z1002, Promega, Madison, WI, USA) according to the manufacturer’s instructions. Brief deparaffinization was performed using mineral oil at 80°C for 1 min, followed by the addition of 100 μL of lysis buffer, centrifugation at 10,000 ×g for 15 s, addition of 10 μL of proteinase K, incubation at 56°C for 15 min, and incubation at 80°C for 1 h. Then, 30 μL of DNase mix was added directly to the lower phase and incubated for 30 min at room temperature. In total, 325 μL of BL buffer and 200 μL of isopropanol (100%) were added to the lysed sample before vortexing and centrifugation at 10,000 ×g for 15 s. The entire lower (aqueous) phase was transferred to a binding column placed in a tube and the assembly was centrifuged at 10,000 ×g for 30 s, washed twice with 500 μL of 1× wash solution, and centrifuged at 10,000 ×g for 30 s. The centrifuge column was dried at 16,000 ×g for 3 min and the RNA was eluted in 50 μL of nuclease-free water centrifuge by centrifuging at 16,000 ×g for 1 min. RNA concentration and purity were determined using a NanoDrop ND-2000 spectrophotometer (Thermo Fisher, Rochester, NY, USA) considering only samples with >1.8 260/280 ratio. RNA integrity was evaluated using agarose gel electrophoresis. All ameloblastoma and dental follicle RNA samples were pooled into a single sample to obtain 2 μg RNA each for the synthesis of cDNA that was used for the microarray.

The cDNA from dental follicle control was labeled with Alexa 555, and that from CAm was labeled with Alexa 647, followed by mixing and hybridization with the GeneChip Human Mapping 10K Array (CHIP H10KA_07_38, AFFYMETRIX, Santa Clara, CA). The microarray data quantification of the chip images was analyzed using genArise software. A Z-Score cutoff of >2.0 was used to determine the upregulated genes. Then, R analysis was performed considering a Benjamini & Hochberg analysis with a false discovery rate (FDR) p<0.05 as significant. The microarray service of the Microarray Unit of the Cellular Physiology Institute of UNAM was used.¹⁵

Functional enrichment analyses

The target gene list was submitted to DAVID 6.8, available online: https://david.ncifcrf.gov.¹⁶ Differentially expressed genes (DEGs), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and Gene Ontology (GO) enrichment analyses were conducted to analyze the functions of the candidate, mainly including biological process (BP), molecular function (MF), and cellular component (CC). Only elements of statistical significative with a p<0.05 were selected.

Protein-protein interaction (PPI)

The Search Tool for the Retrieval of Interacting Genes (STRING); version 11.0, http://string-db.org/ database was used to predict the protein-protein interaction (PPI) networks of the DEGs.¹⁷ Then, the Cytoscape software was used to analyze the interaction with a combined score of >0.4 (http://cytoscape.org). Finally, the plugin molecular complex detection (MCODE) was used to screen the most significant module in the PPI networks with the MCODE score >, degree cutoff=2, node score cutoff=0.2, k-core 2, and max depth=100.

Selection and analyses of hub genes

For the selection of the hub genes, those clustered with MCODE score ≥ 2.5 were selected, and then the effect of the hub genes on overall survival and disease-free survival was analyzed using the Kaplan-Meier plotter (KM plotter, http://kmplot.com/analysis) by adjusting the follow-up threshold to 60 months.¹⁸ To validate these hub genes, we used the Gene Expression Profiling Interactive Analysis (GEPIA Online: http://gepia.cancer-pku.cn/index.html) website to analyze data pertaining to RNA expression from thousands of samples from the Genotype-Tissue Expression (GTEx) and The Cancer Genome Atlas (TCGA) projects. To mimic the behavior of ameloblastoma to the maximum extent of the candidate genes, both analyses were adjusted for head-neck squamous carcinoma.¹⁹

Real-time quantitative Reverse Transcription PCR (RT-qPCR)

Eight additional conventional ameloblastoma FFPE samples and three dental follicles as control were obtained from the Oral Medicine and Pathology Department and the Histopathological archive of the Oral and Maxillofacial Pathology specialty, ENES Leon, UNAM, their histological pattern was determined by two oral pathologists (Supplementary Table 2). The total RNA was obtained by using the ReliaPre FFPE Total RNA Miniprep System. Quantitative Reverse Transcription-PCR was performed using GoTaq 1-Step RT-qPCR System (Promega, Madison, Wisconsin, USA) according to the manufacturer’s instructions, and the data collection was performed on the ABI PRIS 7000 Sequence Detection Systems (Waltham, Massachusetts, USA). The primers for PDGFA were 5’-TTCTGGCTTTGTGTTTCTCCCTTA -3’ (sense) and 5’-TACGATTGGTTGACGCATAGTTCT-3’ (antisense); and for IL2RA were 5’- CAGGAACAGAAGGATGAATGAG-3’ (sense) and 5’- CCAATTAGTAACGCACAGGTAA-3’ (antisense); GAPDH primers were 5´-ACCACAGTCCATGCCATCAC-3´ (sense) and 5´-TCCACCACCCTGTTGCTGTA-3´ (antisense). Relative expression was computed using the 2^-(∆∆Ct) method.

Immunohistochemistry assay and interpretation

The same eight additional conventional ameloblastoma FFPE samples and three dental follicles were employed for immunohistochemical analysis. Three poly-L-lysine-treated slides with 4-µm sections were obtained from each sample. The slides were deparaffinated and rehydrated conventionally in xylene and alcohol washes. Antigenic retrieval was performed in 10 mM citrate buffer in microwave heat (700 W for 3 min and 30 s). Endogenous peroxidase was inhibited with hydrogen peroxide at 3% for 20 min (Sigma-Aldrich, St. Louis, MO, USA). The background was blocked with 100 μl of 2% albumin for 20 min (Sigma-Aldrich, St. Louis, MO, USA), and 100 μl of 0.2% X-100 Triton was then added for 20 min (Sigma-Aldrich, St. Louis, MO, USA). The slides were incubated with primary antibodies for PDGFA (Santa Cruz Biotechnology, sc-9974, Santa Cruz, California, USA) and IL2RA (Santa Cruz Biotechnology, sc-665) at adjusted concentrations of 1:200 overnight at 4°C. Negatives controls were produced by the omission of primary antibodies and substituted with PBS as previously reported.¹⁴

After incubation, primary antibodies were removed and Immunodetector Biotinylated Link was subsequently added and incubated for 10 min (Bio SB, BSB 0007, Santa Barbara, California, USA). The slides were then incubated with Immunodetector HRP label for 10 min. Immunocomplexes were visualized via diaminobenzidine (DAB) incubation for 1 min, and the slides were counterstained for 1 min with Mayer’s hematoxylin. The slides were observed using a Leica DM750 microscope.

For the analysis of each marker, photomicrographs of 5 fields were obtained at 400× magnification from each sample using a Leica ICC50 HD camera. The intensity of staining (optical density) was obtained using the ImageJ software (NIH, Bethesda MD, USA); calibrating the quantification to establish the scale of optical density at: 0-0.9/negative, 1-1.9/mild, 2-2.9/moderate, and >3/intense.

Statistical analysis

Statistical software of SPSS 26.0 (SPSS Inc., Chicago, IL, USA) was employed for the statistical analysis of data. The clinical-demographic data as age, gender, anatomic zone, and histological pattern were analyzed for descriptive distribution and central tendency (mean±standard deviation; x±sd). For the 2-ΔΔCt method, first we estimate the average Ct values for any technical replicates, then we estimate the delta Ct for each sample by the formula:

$\Delta Ct=Ct\text{(gene of interest)}-Ct\text{(housekeeping gene)}$

As a calibrator, we estimate average Ct. To calculate the ∆∆Ct values we employed the following formula:

$\Delta \Delta Ct=\Delta Ct\text{(Sample)}-\Delta Ct\text{(calibrator average)}$

The next step was doing the 2^-(∆∆Ct) of each sample to obtain the average 2^-(∆∆Ct) of control and ameloblastoma. The Shapiro-Wilk test was performed for immunohistochemistry results to determinate the normal distribution of data obtaining a p>0.05 not rejecting the null hypostasis, then an independent sample t-test was employed for the comparison between groups, considering p<0.05 as statistically significant.

Results

The gender distribution was nine males and six females. The mean age was 37.8±17.7 years old. In total, 12 conventional ameloblastoma presented follicular and three plexiform patterns. All specimens were located in the mandible (Supplementary Table 1).

Although the quality of the RNA samples was high, the amount of RNA was insufficient to perform independent microarray analyses. Thus, the RNA samples from ameloblastoma and dental follicles were pooled to obtain ameloblastoma and dental follicle control groups, respectively.

Identification of DEG gene ontology and KEGG pathway analysis

In total, 376 upregulated genes were identified. All genes were analyzed by the DAVID enrichment software, and the results of GO analysis indicated that 1) biological processes (BP) were particularly enriched in positive regulation of cell division, regulation of phosphatidylinositol 3-kinase signaling, branching involved in salivary gland morphogenesis, positive regulation of MAP kinase activity, positive regulation of mesenchymal cell proliferation, negative regulation of transcription from RNA polymerase II promoter, regulation of branching involved in salivary gland morphogenesis by mesenchymal-epithelial signaling, response to wounding, angiogenesis, fibroblast growth factor receptor signaling pathway, positive regulation of peptidyl-tyrosine phosphorylation, immune response, epidermis development, lung-associated mesenchyme development, Notch signaling pathway, response to drugs, induction of positive chemotaxis, epithelial tube branching involved in lung morphogenesis, negative regulation of apoptotic signaling pathway, positive regulation of GTPase activity, peptidyl-tyrosine phosphorylation, and protein localization to cell surface; 2) molecular functions (MF) were enriched in phosphatidylinositol-4,5-bisphosphate 3-kinase activity, Ras guanyl-nucleotide exchange factor activity, growth factor activity, protein tyrosine kinase activity, protein heterodimerization activity, heparin binding, protein phosphatase binding, protein domain specific binding, transcription factor binding, alpha-actinin binding, protein binding, platelet-derived growth factor receptor binding, enzyme binding, protein dimerization activity, and TATA-binding protein (TBP)-class protein binding; 3) cell components (CC) were enriched in the extracellular region, membrane, extracellular space, transcription factor complex, and flotillin complex (Table 1).

Table 1
GO analysis of differential expressed genes associated with ameloblastoma

Ontology	Term	Count	P-Value	FDR
BP	GO:0051781~positive regulation of cell division	6	0.00000412	0.0022802
BP	GO:0014066~regulation of phosphatidylinositol 3-kinase signaling	22	0.0000313	0.00865839
BP	GO:0060445~branching involved in salivary gland morphogenesis	4	0.000285	0.03155322
BP	GO:0043406~positive regulation of MAP kinase activity	15	0.000344	0.03171901
BP	GO:0002053~positive regulation of mesenchymal cell proliferation	4	0.00164912	0.13028075
BP	GO:0000122~negative regulation of transcription from RNA polymerase II promoter	8	0.00592798	0.36424165
BP	GO:0009611~response to wounding	4	0.00938832	0.45943785
BP	GO:0001525~angiogenesis	5	0.0148601	0.5106253
BP	GO:0008543~fibroblast growth factor receptor signaling pathway	4	0.01552901	0.5106253
BP	GO:0050731~positive regulation of peptidyl-tyrosine phosphorylation	4	0.01552901	0.5106253
BP	GO:0006955~immune response	6	0.0160365	0.5106253
BP	GO:0008544~epidermis development	4	0.01662072	0.5106253
BP	GO:0007219~Notch signaling pathway	4	0.02921914	0.80790927
BP	GO:0042493~response to drug	5	0.03324101	0.8244835
BP	GO:0043547~positive regulation of GTPase activity	6	0.04120035	0.9113518
BP	GO:0018108~peptidyl-tyrosine phosphorylation	4	0.04906818	0.96887935
CC	GO:0005576~extracellular region	17	0.000648	0.04653223
CC	GO:0016020~membrane	19	0.000745	0.04653223
CC	GO:0005615~extracellular space	15	0.00258885	0.10786879
CC	GO:0005667~transcription factor complex	9	0.00927696	0.28990486
CC	GO:0016600~flotillin complex	7	0.02007113	0.50177824
MF	GO:0046934~phosphatidylinositol-4.5-bisphosphate 3-kinase activity	7	0.0000111	0.00162983
MF	GO:0005088~Ras guanyl-nucleotide exchange factor activity	7	0.000126	0.00927298
MF	GO:0008083~growth factor activity	7	0.000468	0.02292204
MF	GO:0004713~protein tyrosine kinase activity	6	0.00329741	0.1082232
MF	GO:0046982~protein heterodimerization activity	8	0.00368106	0.1082232
MF	GO:0008201~heparin binding	6	0.00552787	0.13543269
MF	GO:0019903~protein phosphatase binding	3	0.00883842	0.18560683
MF	GO:0019904~protein domain specific binding	6	0.01132918	0.2081737
MF	GO:0008134~transcription factor binding	6	0.02576728	0.40442461
MF	GO:0051393~alpha-actinin binding	4	0.0288818	0.40442461
MF	GO:0005515~protein binding	29	0.03026307	0.40442461
MF	GO:0005161~platelet-derived growth factor receptor binding	4	0.03325244	0.40734242
MF	GO:0019899~enzyme binding	6	0.03860665	0.43655217
MF	GO:0046983~protein dimerization activity	5	0.0447642	0.45324775
MF	GO:0017025~TBP-class protein binding	4	0.04624977	0.45324775

KEGG analysis data (Table 2) showed that the upregulated genes were enriched in pathways in cancer, PI3K-Akt signaling, and Jak-STAT signaling pathways.

Table 2
KEGG pathway analysis of differential expressed genes associated with ameloblastoma

Term	Count	P-Value	Genes	FDR
hsa05215:Prostate cancer	7	0.00337144	PDGFA, CTNNB1, ERBB2, FGFR1, FOS, DVL2, FZD4	0.15462207
hsa04015:Rap1 signaling pathway	6	0.00554887	PDGFA, CTNNB1, VEGFD, FGF7, FGFR1, PAK1	0.15462207
hsa04151:PI3K-Akt signaling pathway	9	0.005947	IL2RA, GH1, PDGFA, VEGFD, FGF7, FGFR1, NDP, IKBKG, AKT3	0.15462207
hsa05200:Pathways in cancer	12	0.01022767	PDGFA, CTNNB1, VEGFD, FGF7, FGFR1, ERBB2, POU5F1, GREM1, NBL, FZD4, GRB2	0.19943962
hsa04630:Jak-STAT signaling pathway	9	0.01342465	IL2RA, GH1, IL13, IL11RA, MPL, THPO, THBD, STAT1, JAK3	0.2094245
hsa05218:Melanoma	4	0.02504124	PDGFA, FGF7, FGFR1, PIN1	0.27903096
hsa04520:Adherens junction	3	0.02504124	CTNNB1, ERBB2, CD70	0.27903096
hsa04510:Focal adhesion	8	0.03371729	PDGFA, CTNNB1, FGFR1, VEGFD, ERBB2, THPO, THBD, MRC2	0.3287436
hsa04014:Ras signaling pathway	5	0.04260335	PDGFA, VEGFD, FGF7, FGFR1, TIAM	0.36922906

PPI network and modular analysis

A total of 31 DEGs, including 29 nodes and 68 edges, were imported into the PPI network complex. We then applied Cytoscape MCODE for further analysis, and the results showed 25 nodes and 68 edges, with 20 clustered genes with an MCODE score >2.5 (Figure 1).

Figure 1
PPI of the differentially expressed genes constructed using STRING online database and Cytoscape software analysis. a) A total of 31 differentially expressed genes were identified in the network. b) Module analysis via Cytoscape software (degree cutoff=2, node score cutoff=0.2, k-core=2, and max depth=100)

Analysis of core genes using Kaplan-Meier plotter and GEPIA

The Kaplan-Meier plotter was used to identify the survival data for these 20 clustered genes. Only nine genes were significantly associated with poor survival (Figure 2). GEPIA was used to validate these nine genes and led to the identification of two genes (platelet-derived growth factor A (PDGFA) and interleukin 2 receptor subunit alpha (IL2RA) with significant correlation (Figure 3).

Figure 2
A) Prognostic information of the 20 core genes. Kaplan-Meier plotter online tool was used to analyze the prognostic information and nine genes were found to be significantly associated with survival rate (*p<0.05). B) Validation of the significant genes by GEPIA. The significant genes expressed in patients with ameloblastoma were compared to those in healthy individuals. Only platelet-derived growth factor A (PDGFA) and interleukin 2 receptor subunit alpha (IL2RA) showed significant differential expression (*p<0.05)

Figure 3
Selection of candidate genes

IL2RA and PDGFA gene expression and immunohistochemistry analysis

Three samples showed a follicular pattern and five were plexiform (Supplementary Table 2). The gene expressions of IL2RA and PDGFA in conventional ameloblastoma were higher than those in the dental follicle in the 2^-(∆∆Ct) method relative quantification with 362±66 and 419±33 measure units, respectively. The immunoexpression analysis showed that IL2RA presented an intense expression in the parenchyma and stroma of CAm, especially in the follicular pattern. The PDGFA showed a moderate to mild immunoexpression in plexiform and follicular patterns respectively, predominantly in the parenchyma, however, there was no significant difference related to histological pattern (Figure 4).

Figure 4
Immunoexpression analysis. a) and b) showed IL2RA intense immunoexpression in stroma and parenchyma of conventional ameloblastoma. c) and d) showed more selective immunoexpression of PDGFA for parenchyma

Discussion

Conventional ameloblastoma is a benign epithelial odontogenic tumor that is frequently diagnosed in young adults with a median age of 35 years without any gender-specific trend. CAm often progresses slowly but is locally invasive. Untreated tumors resorb the cortical plate bone and extend into the adjacent tissue.⁵,²⁰ Our samples used for microarray were obtained from 15 patients (six females and nine males) with a mean age of 37.2±17.8 years, which is consistent with that reported previously.²⁰ However, it has been proposed that when mutation BRAF V600E is present, the presentation age is earlier than the wildtype genotype.² In Mexico, CAm is commonly diagnosed in advanced stages due to the absence of symptoms and low prevalence, which results in detrimental effects on the bone as described above, thus complicating the patient’s treatment and prognosis.³ Increasing our knowledge of the mechanisms underlying the pathogenesis of ameloblastoma is necessary to improve the prognosis and treatment approaches of patients with CAm, as well as the number of possible therapeutic targets. Signaling pathways of WNT, Akt, and FGFR1, effects on bone remodeling by RANK-RANKL and OPG, degradation of extracellular matrix by MMP, and mutations in BRAF and SMO are all molecular events associated with CAm pathogenesis, however, the high-throughput assay could improve this knowledge.²,⁷,⁸

The application of bioinformatics methods on microarray profile datasets is an important strategy to identify more useful therapeutic and/or prognostic biomarkers of ameloblastoma. A total of 376 upregulated genes (log FC>2) were identified. We performed GO function and KEGG pathway enrichment analysis using the DAVID platform to understand the functional relevance of these DEGs and found that, for biological processes, molecular functions and cell components of the upregulated genes were particularly enriched in positive regulation of cell division, positive regulation of MAP kinase activity, positive regulation of mesenchymal cell proliferation, growth factor activity, protein tyrosine kinase activity, and other elements related to cancer development. These data reinforce that although ameloblastoma is classified as a benign tumor, molecularly it emulates the patterns observed in malignant neoplasms, such as gastric cancer, nasopharyngeal carcinoma, breast, or prostate cancer.²¹-²⁴Hu, et al.¹⁰ (2016) in their pathway and gene enrichment analysis showed that inflammation, MAP kinase, and cell cycle regulation are differentially expressed, however when separate the expression of pre-secretory ameloblast and odontoblast contrasting differences was observed. In high-throughput analysis as microarrays, the separate analysis of tumoral parenchyma and stroma via laser capture microdissection could bring information that is validated when immunohistochemistry is performed for protein determination and Nanostring gene expression analysis.¹⁰,¹¹ In our approach, the RT-qPCR and immunohistochemistry were performed, despite not having separated the parenchyma from the stroma since the microarray analysis. Our methodology provides us the possibility of observing and analyzing the tumor microenvironment in a comprehensive way, in which the most relevant elements stand out, and after their validation with a more sensitive assay, as above mentioned, the results could be correlated directly with clinical variables.

We built a PPI network complex with 29 nodes and 68 edges using the STRING online database and Cytoscape software. Twenty upregulated hub genes were screened from the PPI network complex by Cytotype MCODE analysis. Moreover, the Kaplan-Meier plotter analysis revealed that nine of the 20 genes were significantly associated with poor survival. Validation of these nine genes led to the identification of two genes that showed high and significative expression (P<0.05) compared to normal samples by GEPIA analysis (PDGFA and IL2RA) and these represent potentially new effective targets to improve the prognosis or treatment of CAm. Our bioinformatic protocol or algorithm is based on the genome-proteome-clinical utility premise; using GO and KEGG to assess the enriched genome panorama; with protein-protein interaction to estimate the proteome status; and with Kaplan Meier and GEPiA analysis to search the relationship with survival. The hub genes were validated by RT-qPCR and immunohistochemistry to confirm the high presence and observe the distribution in tumoral parenchyma and stroma. Although these techniques only allow the detection of a single candidate per test, they offer greater sensitivity; and the possibility to establish their association with cytological or histological criteria.

PDGF promotes cell proliferation, survival, and migration. Alterations in signaling have been observed in cancer, fibrosis, and atherosclerosis. PDGF is an important factor in ameloblastoma pathogenesis and the expression of the PDGF chain is higher in ameloblastic tumors than in tooth germs, and alongside its cognate receptor (PDGFRA) it is expressed at a variable level in ameloblastomas.²⁵,²⁶

PDGF signaling is important for the growth and differentiation of stem cells, particularly mesenchymal cells. Their dimerization promotes autophosphorylation at ten sites, which can interact with SH2-domain-containing signaling proteins. Activated signaling proteins include phospholipase C, PI3K, Grb2, and others. Binding these proteins leads to the activation of several signaling pathways as MAP kinase pathways, PI3-kinase-Akt, and PLC pathways.²³If we consider PDGFA as a possible therapeutic target, there are three main approaches to inhibit the PDGF/PDGFR pathway: 1) sequestering the ligand with neutralizing antibodies, 2) blocking the receptor with receptor-specific antibodies or small molecules inhibitors, and 3) blocking the kinase activity of PDGFR using low molecular weight inhibitors. Olaratumab, nilotinib, dasatinib, ponatinib, sunitinib, imatinib, and other anti-PDGF drugs have shown significant clinical results in many malignancies. A similar strategy has been developed for BRAF V600E mutation-positive ameloblastomas, in which treatment with vemurafenib showed significant clinical results.²⁷ Many reports postulate that the aggressiveness and frequency of positive BRAF V600E ameloblastoma are high.²⁸-³¹

Cytokines are major mediators of the immune response. Cytokine ligands and receptors control various cellular functions, including proliferation, differentiation, and cell survival/apoptosis of leukocytes; however, they are also involved in many pathophysiological processes. The interleukin-2 receptor is involved in the regulation of immune tolerance by controlling regulatory T cell (TREG) activity. The interleukin 2 (IL2) receptor alpha (IL2RA), beta (IL2RB), and gamma chain (IL2RG) constitute high-affinity IL2 receptors. Homodimeric alpha chains (IL2RA) result in a low-affinity receptor, whereas homodimeric beta (IL2RB) chains produce a medium-affinity receptor.³² Ameloblastoma cells and surrounding stromal cells, such as fibroblasts, may contribute to ameloblastoma pathogenesis.⁸Chantravekin and Koontongkaew conducted a co-culture test and three-dimensional organotypic culture and showed the role of the stroma on the tumor parenchyma, as fibroblasts associated with ameloblastomas modulated tumor development, promoting the proliferation and invasion of tumor cells.³³ It has been hypothesized that ameloblastoma cells and stromal fibroblasts may be reciprocally activated via cytokines, such as IL-6, IL-8, IL-1a, and recently IL33, to create a tumoral microenvironment that promotes tumor formation.⁸,³⁴ Damoiseaux³⁵ (2020) has reported that the IL2RA fraction can function in a diverse way to lead to leukocyte activation via paracellular or even soluble forms, which can affect the functionality of cells such as cytotoxic CD4 and CD8 T lymphocytes at the tumor level. This is the first report to identify IL2RA as a possible participant in the mechanism underlying the development of ameloblastoma. When we analyzed the KM plot result, we observed that patients with a high level of IL2RA have a greater probability of survival. If we correlate this result with immunoexpression results, immunomodulation by IL2RA is present, but additional studies would need to verify the mechanism that conducts it, as well as the result of PDGFA, since due to its greater tendency to express itself in tumor parenchyma, it becomes a direct target. Nevertheless, the validation analysis of IL2RA and PDGFA by immunohistochemistry reinforces the concept that the parenchyma and stroma relationship is a necessary feature that must be considered to improve our understanding, in order to develop better therapeutic strategies. Zhang, et al.³⁶ (2022) suggested in a bioinformatic analysis that macrophages could infiltrate the ameloblastoma and participate in their pathogenesis.³⁶ That could be a relationship mainly with IL2RA, however, to prove this immunological relationship their validation is necessary.

Taken together, our bioinformatic analysis identified two hub genes (PDGFA and IL2RA) between CAm and normal dental follicles. The results suggested that these genes play key roles in the pathogenesis, progression, and prognosis of CAm. The main limitation to postulating PDGFA and IL2RA as therapeutic targets is the verification of their reach in cellular or animal models, in which the biological behavior could be measured and correlated. For this reason, identifying how we can affect CAm in these specific targets to provide useful information on these new biomarkers is necessary.

Acknowledgments

This study was supported with institutional resources from PAPIIT IN226720, National Autonomous University of Mexico.

Reference

• ¹
  - Zheng L, Ehardt L, Mc Alpm B, About I, Kim D, Papagerakis S, et al. The tick tock of Odontogenesis. Exp Cell Res. 2014;325(2):83-9. doi: 10.1016/j.yexcr.2014.02.007
• ²
  - EI-Naggar AK, Chan JKC, Grandis JR, Takata T, Slootweg PJ, editors. WHO Classification of Head and Neck Tumours. 4th ed. Lyon: International Agency for Research on Cancer; 2017. 296 p.
• ³
  - Gaitán-Cepeda LA, Quezada-Rivera D, Tenorio-Rocha F, Leyva-Huerta ER. Reclassification of odontogenic keratocyst as tumour: impact on the odontogenic tumours prevalence. Oral Dis. 2010;16(2):185-7. doi: 10.1111/j.1601-0825.2009.01625.x
• ⁴
  - Gomes CC, Duarte AP, Diniz MG, Gomez RS. Review article: current concepts of ameloblastoma pathogenesis. J Oral Pathol Med. 2010;39(8):585-91. doi: 10.1111/j.1600-0714.2010.00908.x
• ⁵
  - Laborde A, Nicot R, Wojcik T, Ferri J, Raoul G. Ameloblastoma of the jaws: management and recurrence rate. Eur Ann Otorhinolaryngol Head Neck Dis. 2017;134(1):7-11. doi: 10.1016/j.anorl.2016.09.004
• ⁶
  - Bell DW. Our changing view of the genomic landscape of cancer. J Pathol. 2010;220(2):231-43. doi: 10.1002/path.2645
• ⁷
  - Brown NA, Betz BL. Ameloblastoma: a review of recent molecular pathogenetic discoveries. Biomark Cancer. 2015;7(Suppl 2):19-24. doi: 10.4137/BIC.S29329
• ⁸
  - Fuchigami T, Ono Y, Kishida S, Nakamura N. Molecular biological findings of ameloblastoma. Jpn Dent Sci Rev. 2021;57:27-32. doi: 10.1016/j.jdsr.2020.12.003
• ⁹
  - Heikinheimo K, Jee KJ, Niini T, Aalto Y, Happonen RP, Leivo I, et al. Gene expression profiling of ameloblastoma and human tooth germ by means of a cDNA microarray. J Dent Res. 2002;81(8):525-30. doi: 10.1177/154405910208100805
• ¹⁰
  - Hu S, Parker J, Divaris K, Padilla R, Murrah V, Wright JT. Ameloblastoma phenotypes reflected in distinct transcriptome profiles. Sci Rep. 2016;6:30867. doi: 10.1038/srep30867
• ¹¹
  - DeVilliers P, Suggs C, Simmons D, Murrah V, Wright JT. Microgenomics of ameloblastoma. J Dent Res. 2011;90(4):463-9. doi: 10.1177/0022034510391791
• ¹²
  - Ramírez-Martínez CM, Jacinto-Alemán LF, Cruz-Hervert LP, Portilla-Robertson J, Leyva-Huerta ER. Bioinformatic analysis for mucoepidermoid and adenoid cystic carcinoma of therapeutic targets. Vaccines (Basel). 2022;10(9):1557. doi: 10.3390/vaccines10091557.
• ¹³
  - Universidad Nacional Autónoma de México. Facultad de Odontología. Aviso de privacidad integral [Internet]. Ciudad de México: UNAM; 2019 [cited 2021 Feb 28]. Available from: http://www.odonto.unam.mx/sites/default/files/inline-files/1.AP_Integral-Fac_Esc_Inst_Centros%20%20_0.pdf
  » http://www.odonto.unam.mx/sites/default/files/inline-files/1.AP_Integral-Fac_Esc_Inst_Centros%20%20_0.pdf
• ¹⁴
  - Trejo-Remigio DA, Jacinto-Alemán LF, Leyva-Huerta ER, Navarro-Bustos BR, Portilla-Robertson J. Ectodermal and ectomesenchymal marker expression in primary cell lines of complex and compound odontomas: a pilot study. Minerva Stomatol. 2019;68(3):132-41. doi: 10.23736/S0026-4970.19.04166-9
• ¹⁵
  - Tecalco-Cruz AC, Macías-Silva M, Ramírez-Jarquín JO, Méndez-Ambrosio B. Identification of genes modulated by interferon gamma in breast cancer cells. Biochem Biophys Rep. 2021;27:101053. doi: 10.1016/j.bbrep.2021.101053
• ¹⁶
  - National Cancer Institute. Frederick National Laboratory for Cancer Research. Laboratory of Human Retrovirology and Immunoinformatics. DAVID Bioinformatics Resources 6.8 [Internet]. Frederick, Laboratory of Human Retrovirology and Immunoinformatics; c2021 [cited 2021 Sept 17]. Available from: https://david.ncifcrf.gov/
  » https://david.ncifcrf.gov/
• ¹⁷
  - String Consortium. STRING 11.0 [Internet]. Cambridgeshire: ELIXIR; c2021 [cited 2021 Sept 17]. Available from: https://string-db.org/
  » https://string-db.org/
• ¹⁸
  - Nagy A, Munkacsy G, Gyorffy B. Pancancer survival analysis of cancer hallmark genes. Sci Rep. 2021;11(1):6047. doi: 10.1038/s41598-021-84787-5
• ¹⁹
  - Tang Z, Li C, Kang B. GEPIA - Gene Expression Profiling Interactive Analysis [Internet]. Beijing: Peking University, Zhang Lab; c2017 [cited 2021 Sept 17]. Available from: http://gepia.cancer-pku.cn/index.html
  » http://gepia.cancer-pku.cn/index.html
• ²⁰
  - Masthan KM, Anitha N, Krupaa J, Manikkam S. Ameloblastoma. J Pharm Bioallied Sci. 2015;7(Suppl 1):S167-70. doi: 10.4103/0975-7406.155891
• ²¹
  - Ding C, Zhang Q, Jiang X, Wei D, Xu S, Li Q, et al. The analysis of potential diagnostic and therapeutic targets for the occurrence and development of gastric cancer based on bioinformatics. Comput Math Methods Med. 2022:4321466. doi: 10.1155/2022/4321466
• ²²
  - Song Y, Feng T, Cao W, Yu H, Zhang Z. Identification of key genes in nasopharyngeal carcinoma based on bioinformatics analysis. Comput Intell Neurosci. 2022:9022700. doi: 10.1155/2022/9022700
• ²³
  - Mehraj U, Alshehri B, Khan AA, Bhat AA, Bagga P, Wani NA, et al. Expression pattern and prognostic significance of chemokines in breast cancer: an integrated bioinformatics analysis. Clin Breast Cancer. 2022;22(6):567-78. doi: 10.1016/j.clbc.2022.04.008
• ²⁴
  - Khan MM, Mohsen MT, Malik MZ, Bagabir SA, Alkhanani MF, Haque S, et al. Identification of potential key genes in prostate cancer with gene expression, pivotal pathways and regulatory networks analysis using integrated bioinformatics methods. Genes (Basel). 2022;13(4):655. doi: 10.3390/genes13040655
• ²⁵
  - Kumamoto H, Ooya K. Immunohistochemical detection of insulin-like growth factors, platelet-derived growth factor, and their receptors in ameloblastic tumors. J Oral Pathol Med. 2007;36(4):198-206. doi: 10.1111/j.1600-0714.2007.00516.x
• ²⁶
  - Sulzbacher I, Wick N, Pichlhofer B, Mazal PR. Expression of platelet-derived growth factor-AA and platelet-derived growth factor-alpha receptor in ameloblastomas. J Oral Pathol Med. 2008;37(4):235-40. doi: 10.1111/j.1600-0714.2008.00637.x
• ²⁷
  - Papadopoulos N, Lennartsson J. The PDGF/PDGFR pathway as a drug target. Mol Aspects Med. 2018;62:75-88. doi: 10.1016/j.mam.2017.11.007
• ²⁸
  - Fernandes GS, Girardi DM, Bernardes JPG, Fonseca FP, Fregnani ER. Clinical benefit and radiological response with BRAF inhibitor in a patient with recurrent ameloblastoma harboring V600E mutation. BMC Cancer. 2018;18(1):887. doi: 10.1186/s12885-018-4802-y
• ²⁹
  - Seki-Soda M, Sano T, Ito K, Yokoo S, Oyama T. An immunohistochemical and genetic study of BRAFV600E mutation in Japanese patients with ameloblastoma. Pathol Int. 2020;70(4):224-30. doi: 10.1111/pin.12899
• ³⁰
  - Oh KY, Cho SD, Yoon HJ, Lee JI, Ahn SH, Hong SD. High prevalence of BRAF V600E mutations in Korean patients with ameloblastoma: clinicopathological significance and correlation with epithelial-mesenchymal transition. J Oral Pathol Med. 2019;48(5):413-20. doi: 10.1111/jop.12851
• ³¹
  - Derakhshan S, Aminishakib P, Karimi A, Saffar H, Abdollahi A, Mohammadpour H, et al. High frequency of BRAF V600E mutation in Iranian population ameloblastomas. Med Oral Patol Oral Cir Bucal. 2020;25(4):e502-7. doi: 10.4317/medoral.23519
• ³²
  - Ohue Y, Nishikawa H. Regulatory T (Treg) cells in cancer: can Treg cells be a new therapeutic target? Cancer Sci. 2019;110(7):2080-9. doi: 10.1111/cas.14069
• ³³
  - Chantravekin Y, Koontongkaew S. Effects of ameloblastoma-associated fibroblasts on the proliferation and invasion of tumor cells. J Cancer Res Ther. 2014;10(4):1082-7. doi: 10.4103/0973-1482.138005
• ³⁴
  - Chen L, Wang G, Qiao X, Wang X, Liu J, Niu X, et al. Downregulated miR-524-5p Participates in the tumor microenvironment of ameloblastoma by targeting the interleukin-33 (IL-33)/suppression of tumorigenicity 2 (ST2) axis. Med Sci Monit. 2020;26:e921863. doi: 10.12659/MSM.921863
• ³⁵
  - Damoiseaux J. The IL-2 - IL-2 receptor pathway in health and disease: the role of the soluble IL-2 receptor. Clin Immunol. 2020;218:108515. doi: 10.1016/j.clim.2020.108515
• ³⁶
  - Zhang Z, Peng Y, Dang J, Liu X, Zhu D, Zhang Y, et al. Identification of key biomarkers related to epithelial-mesenchymal transition and immune infiltration in ameloblastoma using integrated bioinformatics analysis. Oral Dis. Forthcoming 2022. doi: 10.1111/odi.14173

• Data availability statement
  The datasets generated or analyzed during this study are available in the SciELO Data repository, [DOI 10.48331/scielodata.Z2S8X9].

Publication Dates

• Publication in this collection
  23 Jan 2023
• Date of issue
  2022

History

• Received
  9 Aug 2022
• Reviewed
  18 Nov 2022
• Accepted
  21 Nov 2022