Identification of differentially expressed genes (DEGs)
A total of 65 human GIST samples without preoperative imatinib treatment were included in dataset GSE13675. A cohort of 47 samples with KIT mutations, consisting of 16 advanced (high-risk) GIST patients and 31 non-advanced (low-risk and very low-risk) GIST patients, was used to screen the DEGs using GEO2R. A total of 606 genes (244 upregulated and 362 downregulated) were identified using the cut-off criteria of adj. p < 0.05 and Log2FoldChange > 1. The top 50 up- and down-regulated genes were presented in a heatmap (Figure 1A, supplementary file 1).
Functional annotation of DEGs
To evaluate the biological clustering of DEGs, GO and KEGG analyses for the up- and downregulated DEGs were performed using DAVID. Based on GO analysis, upregulated DEGs were found to be significantly enriched in cell division, sister chromatid cohesion, mitotic nuclear division, microtubule-based movement, and mitotic metaphase plate congression. Cellular component (CC) of the upregulated DEGs was significantly enriched in midbody, kinesin complex, spindle, spindle microtubule, and condensed chromosome kinetochore. Molecular function (MF) of the upregulated DEGs was significantly enriched in microtubule motor activity and microtubule-binding. KEGG analysis showed that the upregulated DEGs were mainly involved in cell cycle and oocyte meiosis. GO analysis revealed that the down-regulated DEGs were significantly enriched in the interferon-gamma-mediated signaling pathway, type I interferon signaling pathway, and antigen processing and presentation of antigens. The CC of the down-regulated DEGs was significantly enriched in the integral component of the lumenal side of the endoplasmic reticulum membrane and ER to Golgi transport vesicle membrane while the MF of the down-regulated DEGs was significantly enriched in peptide antigen binding and MHC class II receptor activity. KEGG analysis further revealed that the down-regulated DEGs were mainly involved in Graft-versus-host disease, Type I diabetes mellitus, Allograft rejection, Antigen processing and presentation, and Viral myocarditis. These results are presented in Figure 1B (supplementary file 2).
Module analysis through the PPI network of DEGs
To clarify the DEGs functionally, STRING was used to conduct a PPI network, which was composed of 603 nodes and 1768 edges. The PPI enrichment p-value was < 1.0E-16. The PPI network was visualized by Cytoscape and further analyzed by plug-in MCODE. Two most significant modules (module 1 and module 2) were identified and analyzed using GO, KEGG and REACTOME annotations to infer their biological functions. The module 1 (MCODE score = 36.667) was mainly involved in in cell cycle-related biological processes and signaling pathways, while the module 2 (MCODE score = 18.759) was mainly involved in immunological processes and signaling pathways. These results are presented in Figure 2A (supplementary file 2).
Identification of hub genes
To reveal the crucial genes underlying the regulation of GIST progression, we filtered hub genes among DEGs using the plugin-in cytoHubba of Cytoscape. Two algorithms, degree and bottleneck, were applied to weight the DEGs. The degree algorithm calculates the relevance and abundance of genes, while the bottleneck algorithm evaluates key genes positions in an entire regulatory network. According to the degree algorithm, the top 15 hub genes were CDK1, KIF11, KIF2C, CENPE, KIF20A, BUB1, CCNA2, CCNB1, AURKA, MAD2L1, CDCA8, KIF4A, CENPF, NDC80, and KIF23, with their scores ranging from 58 to 48. According to the bottleneck algorithm, the top 15 hub genes were AURKA, FN1, CD44, VEGFA, IL6, HLA-DQB1, HLA-DPA1, CXCL8, NT5E, ANK2, FOXM1, CHEK1, STAT1, CDC25A, and IFIH1, with their scores ranging from 64 to 8. A Venn diagram was used to identify the key hub genes by taking intersection of the two hub gene cohorts. The result showed that AURKA was the only overlapping hub gene (Figure 2B).
Gene Set Enrichment analysis (GSEA)
A total of 15 human gastric GIST samples in GSE47911 includes 6 high-risk cases, 1 intermediate-risk case, 3 low-risk cases, and 5 very low-risk cases. To further verify significant biological processes associated with AURKA expression, GSE47911 gene profiles were divided into two groups after which GSEA was performed based on AURKA expression level. Samples with the highest (25%, 4 samples) and lowest (25%, 4 samples) expression levels were selected for further analysis using GSEA. Cell cycle-related gene sets were associated with elevated AURKA expression (Figure 2C, supplementary file 3).
Correlation between AURKA expression and clinicopathological features in GISTs
To evaluate the clinical significance of AURKA expression in GISTs, AURKA expression levels in 49 GIST tissues were assessed by IHC staining (Figure 3A). The correlations between AURKA expression and clinicopathological features (age, gender, location and risk stratification) were determined (Table 2, supplementary file 4). AURKA expression was closely associated with tumor risk stratification (Figure 3B; P < 0.001). The clinical significance of AURKA expression in GISTs was also evaluated using the data from GSE136755 and the raw data provided by Lagarde et al.[24] (Table 3, Table 4, and Figure 4, supplementary file 5). Findings from the GSE136755 revealed significant associations between AURKA expression and tumor risk stratification (P < 0.001) as well as tumor stage (P < 0.001), while analyses of the raw data from Lagarde et al.[24] also showed a significant association between AURKA expression and tumor risk stratification (P < 0.001) as well as tumor recurrence (P < 0.001) and metastasis (P < 0.001). However, apart from GSE136755 which revealed a significant association between AURKA expression and tumor location (P = 0.018), the data provided by Lagarde et al. did not establish a significant association between AURKA expression and tumor location (P = 0.156).
To determine the prognostic value of AURKA expression in GISTs, Kaplan-Meier survival analysis was performed. The range of observation time was 9 - 79 months. As shown in Figure 3C, GIST patients with elevated AURKA expression exhibited poorer DFS than those with low AURKA expression level (43.25±6.94 months vs 98.48±3.44 months,P < 0.001). COX proportional hazard model showed that AURKA could be used as an independent prognostic marker for GISTs (P = 0.002).
Gene mutation types can predict the responsiveness of GISTs to imatinib. GISTs with KIT exon 11, PDGFRα exon 12 and PDGFRα exon 14 mutations were considered sensitive to imatinib. GISTs with other mutations, such as KIT exon 9, KIT exon 13, KIT exon 14, KIT exon 17, KIT exon 18, PDGFRα exon 18 D842V and KIT/PDGFRα wild type GISTs, were insensitive/resistant to imatinib[1, 3, 30]. In the dataset GSE136755, containing 56 imatinib-sensitive and 7 imatinib-resistant GISTs, there was a weak association between AURKA expression and imatinib-resistant gene mutations (Figure 4A, P = 0.074). Through analysis of the raw data provided by Lagarde et al.[24], which contained 45 imatinib-sensitive and 15 imatinib-resistant samples, it was shown that there was a strong association between AURKA expression and imatinib-resistant gene mutations (Figure 4B, P = 0.018).
AURKA expression patterns in common human malignancies
To determine whether elevated AURKA expression is common in human digestive malignancies, mRNA expression levels of AURKA in stomach carcinoma, liver hepatocellular carcinoma, and colorectal carcinoma were evaluated using data from the GEPIA database. AURKA expression was found to be significantly up-regulated in all the above malignancies by comparison to the normal tissues. Findings from the Oncomine database also supported the upregulated expression of AURKA in most human malignancies (Figure 3D).
AURKA overexpression enhances GIST/T1 cell proliferation, survival, and migration
To assess the biological effects of AURKA expression in GISTs, AURKA was overexpressed in GIST/T1 cells transfected with the AURKA-expressing virus, which was defined as the AURKA overexpression group (AURKA group). Normal GIST/T1 cells (Blank group) and GIST/T1 cells transfected with vacant plasmids (Vector group) were considered as the control groups. The transfection efficiency was determined by the red fluorescence from tdtomato, and quantified by RT-qPCR and western blotting. Figure 5 showed that compared to the Blank and Vector groups, AURKA was overexpressed in GIST/T1 cells transfected with AURKA-expressing virus (the AURKA group). (supplementary file 6)
The CCK8 assay was performed to assess the effect of AURKA overexpression on cell proliferation. Compared to the Blank and Vector groups, the overexpression of AURKA significantly enhanced GIST/T1 cell proliferation (P = 0.018) (Figure 6A). Imatinib treatment significantly inhibited cell proliferation in all the three groups. However, compared to the control groups, AURKA overexpressed cells still showed a relatively higher proliferation rate in the presence of imatinib (P < 0.001, Figure 6A). (supplementary file 6)
We also established that AURKA overexpression markedly suppressed the apoptotic process in GIST/T1 cells (P < 0.001, Figure 6B). Similar result was observed in the circumstance of imatinib administration. Compared to the Blank and Vector groups, AURKA overexpression obviously inhibited cell apoptosis after imatinib administration (P < 0.001, Figure 6B). The results suggested that the AURKA overexpression enhanced resistance of GIST cells to imatinib. (supplementary file 6)
Apart from cell proliferation and apoptosis, the effect of AURKA overexpression on cell migration was also evaluated. The number of AURKA overexpressed cells that passed through the membrane was significantly higher than that of the other two control groups (for AURKA group, rate = 47.00 ± 5.06%; for Blank group, rate = 26.17 ± 3.66%;for Vector group, rate = 30.00 ± 4.15%) in the absence of imatinib (P < 0.001, Figure 6C). However, imatinib significantly suppressed cell migration in all three groups (P < 0.001, Figure 6C) and effectively eliminated the contribution of AURKA overexpression to cell migration (P = 0.169, Figure 6C). (supplementary file 6)