Introduction
Thyroid cancer (TC) is the most common endocrine malignancy, accounting for approximately 90% of all malignancies of the endocrine system (1,2). The morbidity associated with TC is exceedingly high, and an increasing number of asymptomatic TC cases are detected due to the wide use of routine high-resolution ultrasonography (3,4). In the United States, it is estimated that 53,990 cases of TC will be diagnosed in 2018 (2). In fact, TC has become the fifth most common type of cancer affecting women. Despite the fact that TC exhibits an indolent and non-aggressive behavior in the majority of cases, the high incidence rate and several grievous clinical manifestations, such as lymph node metastases, remain unresolved issues (4–6). Based on the histopathology, TC is classified into four main tumor types, including papillary thyroid cancer (PTC), follicular thyroid carcinoma (FTC), medullary thyroid carcinoma (MTC) and anaplastic thyroid cancer (ATC) (7,8). The different subtypes of TC have enormous heterogeneity in terms of their morphological characteristics and prognoses. Although PTC is the predominant type, accounting for 80% of all TCs, it presents a favorable prognostic tendency overall; however, the specific situation of each patient with PTC is not the same (9,10). Therefore, there is an urgent need to further reveal the molecular pathogenesis of TC.
Autophagy is a highly conserved ‘self-devouring’ process that ensures the orderly degradation of the cytoplasmic contents and the recycling of the macromolecular constituents to maintain cellular homeostasis (11). Considering the fundamental roles of autophagy in regulating cellular biological processes, it is no surprise that it plays essential roles in a wide spectrum of physiological conditions and diseases, particularly cancer (12–15). Moreover, exploring the potential mechanisms of action of autophagy would not only expose the mysterious veil of tumorigenesis, but may also aid in the identification of novel therapeutic targets for autophagy in cancer (16). Several studies have reported the probability of harnessing autophagy in the development of cancer (16–18). As regards TC, autophagy-targeted therapy is now considered to be a valuable strategy, as the regulation of the autophagy level can mediate the chemosensitivity and radiosensitivity of TC cells and can markedly affect the regulation of apoptosis (19). For example, Wang et al found that targeting autophagy sensitized BRAF-mutant TC to vemurafenib (20). In addition, it has been suggested that rapamycin can inhibit the invasive ability of TC cells (21). However, these previous studies only examined single oncogenes or suppressor genes in TCs, which reflect the partial functions of autophagy. There are very few systematic comprehensive analyses addressing the clinical significance and potential biological process of autophagy in TC (22,23).
The Human Autophagy Database (HADb) offers an exhaustive and up-to-date list of autophagy-related genes (ARGs) to satisfy the needs of researchers (24). Previously, Zhang et al calculated 74 differentially expressed ARGs, and two of these were associated with the prognosis of patients with glioma (25). However, the clinical values of 234 ARGs in TC remain unclarified. In this study, we described the expression profiles of 234 ARGs, and we identified a class of ARGs with disrupted expression statuses in the onset of TC based on the Cancer Genome Atlas (TCGA) database. Notably, we also identified the biological process pathways enriched by these genes, which may help in the development of more effective treatments. On the basis of the underlying proliferative functions of ARGs, we further emphasized the clinical significance of differentially expressed ARGs in providing novel insight into the clinical management of TC. The present study provides a deeper understanding of the role of autophagy in TC and may help to improve the clinical outcomes of patients with TC.
Materials and methods
Collection of ARGs
Our researchers acquired the ARGs from the HADb. Subsequently, we downloaded the TC gene expression dataset from the TCGA database, from which we extracted the expression levels of the ARGs. TCGA provided 502 TC and 52 non-tumor cases with gene expression profiles. Accordingly, the clinicopathological data of the patients with TC were also downloaded for use in the current study.
Calculation of the differentially expressed ARGs and the functional analysis of their enriched pathways
The edgeR software package was applied in this study to filter and normalize the expression profiles found in the TCGA database, aiming to analyze the differentially expressed ARGs, not only in TC tumors, but also in the adjacent normal tissues. The standard was a fold change of >2 and P<0.05 after correction. We constructed the network of protein-protein interactions (PPIs) with the help of STRING (https://string-db.org/) (26), an online database available to all, and we then input the interacting data derived from STRING to Cytoscape 3.5.1 (27) for visualization. In addition, enrichment analyses were conducted for a better comprehension of the biological functions of the differentially expressed ARGs in TC, which included the analyses of the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Disease Ontology (DO). The enrichment analyses were carried out using the R clusterProfiler package (28) and the results were displayed with the GOplot package (29).
Prognosis index (PI) construction
A univariate Cox analysis was used to determine the genes associated with the survival time of the patients with TC. Subsequently, a univariate Cox analysis was performed to exclude the genes that failed to become independent prognostic biomarkers. Finally, a PI model was constructed with the eligible ARGs based on the formula in which the gene expression was multiplied by the regression coefficient. The Kaplan-Meier (KM) method was used to depict the differences in survival between the high-risk and low-risk groups, and a log-rank analysis was employed for the assessment. In order to investigate the mechanisms of action of ARGs that were consistent with the PI, we obtained the ARG expression, methylation, and copy number variations from the cBioPortal database (30), which provided the data of 508 patients with TC for analysis. We also assessed the ARG splicing events using the TCGA SpliceSeq database (http://bioinformatics.mdanderson.org/TCGASpliceSeq/).
Statistical analysis
SPSS Statistics for Windows version 24.0 (SPSS Inc., Chicago, IL, USA) and R 3.3.1 (https://www.r-project.org/) were utilized for the statistical analyses. R, GraphPad Prism 5 (GraphPad Software Inc., San Diego, CA, USA) and OriginPro 2017 (OriginLab Corp., Northampton, MA, USA) were used for the plotting. Following normalization, all the expression data were transformed to log2 (value +1). We employed an independent t-test to compare the differential expression levels of the ARGs in the TC and adjacent tissues, and to determine the associations between these genes and the clinicopathological characteristics of the patients with TC. Scatter diagrams and boxplots were used to display gene expression. In order to examine the ability of the differentially expressed ARGs to identify tumors, we applied GraphPad software for plotting and calculating the receiver operating characteristic (ROC) curves of each dataset and determining the area under the curve (AUC). The prognostic values of the ARGs were estimated using a univariate regression analysis; subsequently, the prognostic potentials of the ARGs were obtained using a multivariate Cox regression analysis for the construction of the PI model. Correlations between gene mRNA expression levels and methylation degree were examined by Pearson's correlation analysis. A value of P<0.05 was considered to indicate a statistically significant difference.
Results
Differentially expressed ARGs
A total of 234 ARGs were obtained from the HADb database, and 502 TC cases with gene expression profiles were acquired from the TCGA database (Table I). In addition, 52 cases with adjacent normal tissue data were obtained from the TCGA database. After extracting the gene expression levels of 234 cases, we conducted the differential analysis. Conforming to the standards mentioned above, we eventually acquired 31 differentially expressed ARGs, among which 17 were upregulated, while 14 were downregulated (Figs. 1–3). All of the differentially expressed genes are listed in Table II.
Enrichment analyses of the differentially expressed ARGs
Our researchers carried out a series of gene enrichment and pathway analyses, hoping to further investigate the biological functions of autophagy in the development of TC. Using STRING and Cytoscape, we set up the PPI network for these genes (Fig. 4), and discovered the core genes, such as BCL2, JUN and CDKN1A, which cast light on the functions of autophagy in TC and laid the foundation for multi-targeted therapies. The GOplot analysis revealed that in the biological processes, these genes were related to autophagy, as well as to the intrinsic apoptotic signaling pathway, neuron apoptotic process and the response to gamma radiation (Fig. 5A). In terms of the cellular components, these genes participated in the mitochondrial outer membrane, autophagosome and organelle outer membrane functions (Fig. 5B). With regard to the molecular function, these genes played indispensable roles in certain key functions, such as ubiquitin protein ligase binding (Fig. 5C). In addition, KEGG pathway enrichment analysis indicated that these genes were mainly enriched in the pathways relevant to the cell cycle, including autophagy and apoptosis (Fig. 6). Notably, DO analysis (Fig. 7) suggested that these genes were associated with multiple tumors, other than TC, which manifested the fundamental roles of these genes in tumorigenesis and development.
ARG expression in TC and the clinical significance
As the differentially expressed ARGs were largely enriched in tumor-related pathways, it is most likely that these genes could accelerate the development of TC. Therefore, we analyzed the expression of these genes in various clinical parameters and enquired into their associations with the clinical progress. Initially, we verified the expression of these genes in TC. The independent samples t-test and the analysis of the TCGA database indicated that the expression levels of 17 genes were markedly upregulated in the tumor tissues, including CTSB, CDKN1A, BAX, ATG9B, BCL2L1, SPHK1, CX3CL1, DRAM1, DAPK2, ITGA3, BID, ITGB4, DIRAS3, TP63, EVA1A, CDKN2A and SERPINA1. Conversely, 12 genes displayed markedly lower expression levels in the tumor tissues, including ITPR1, JUN, FOS, BCL2, NAJB1, PRKCQ, MAP1LC3C, CCL2, PPP1R15A, GRID2, HSPB8 and TP53INP2 (Fig. 8). Due to the differences between the t-test and edgeR, the IL24 and CXCR4 genes did not show statistical significance. Subsequently, we performed a ROC analysis to determine the ability to differentiate the cancerous and non-cancerous tissues. A total of 28 genes had AUC values >0.7 (Fig. 9), which suggested that these genes were useful in differentiating between cancerous and non-cancerous tissues. Based on the fact that these genes were expressed differentially in thyroid cancer, we further explored the associations between these 31 genes and some of the major clinical parameters: Age (over or under 60 years), sex (male/female), tumor recurrence (yes/no), primary tumor type (multifocal/unifocal), pathological stage (stage III, IV/stage I, II), T stage (stage III, IV/stage I, II), N stage (presence of lymphatic metastasis or not) and M stage (presence of distant metastasis or not). It was revealed that these 31 genes were closely connected with the majority of the parameters, particularly those linked with tumor progression (Table III and Fig. 10). Of these 31 genes, 25 were expressed differentially in the patients with lymphatic metastasis, and 17 displayed marked differences in expression in the T stage. However, we failed to discover the genes closely associated with distant metastasis as these samples were not sufficient for valid research (n=9).
ARG-based PI construction
The results of univariate analysis revealed that 6 ARGs were associated with the prognosis of patients with TC (Figs. 11 and 12). To improve the accuracy and validity of the conclusion, we performed a multivariate Cox regression analysis and ultimately acquired 3 genes: ATG9B, BID and DNAJB1. Accordingly, a PI model was constructed based on these 3 genes: PI= 0.469*ATG9B expression − 0.796*BID expression + 0.782*DNAJB1 expression (Fig. 13). The median PI value was applied to divide the patients with TC into high- or low-risk groups. The KM test was employed in the current study for statistical analysis. The survival analysis revealed that the hazard ratio (HR) of the overall survival (OS) predicted by PI was 4.706 (95% CI, 1.742–12.710; log-rank, P=0.0023), which indicated that the high-risk patients had markedly shorter OS times than their low-risk counterparts (Fig. 14). Additionally, in the high-risk group, the ATG9B and DNAJB1 genes exhibited notably higher expression level s, whereas the BID gene exhibited a significantly lower expression (Fig. 15).
Molecular mechanism of PI
In order to investigate the molecular mechanism of the PI, we analyzed the alternative slicing, methylation, copy number variation and amplification of the ATG9B, BID and DNAJB1 genes. Using the cBioPortal program, we obtained the associations of the expression of these 3 genes with the methylation and copy number variations (Fig. 16). As shown in Fig. 17A, of the 508 TC cases, genetic alterations of these 3 genes were detected in 33 cases. According to the genetic alterations, the patients were categorized into two groups. The KM analysis suggested that the 33 patients with genetic alterations in these 3 genes were likely to have markedly lower survival times than those without the genetic alterations (Fig. 17B, P=0.007). In terms of disease-free survival, despite the fact that there was no significant difference between these two groups, the patients with the genetic alterations were prone to undesirable prognoses (Fig. 17C). Of note, Alternate Terminator was the most significant splicing events of ATG9B (exon 18). Alternate Promoter in exon 2 and exon 3 were the most significant splicing events of BID and DNAJB1 respectively (Fig. 18). However, there was no significant difference between the three alterative splicing events among the different tumor samples.
Discussion
In cellular processes, autophagy is the foundation for homeostatic regulation. An autophagy-perturbed status is common in TC, and it can contribute to tumor progression. The present study highlighted the comprehensive analysis of ARGs in patients with TC. First, we identified the differentially expressed ARGs in TC cases and normal samples. Subsequently, a functional enrichment analysis found that autophagy influenced several tumor-related pathways. We then analyzed the clinical significances of the differentially expressed genes, and found that they were closely related to lymph node metastasis and the tumor stage. It is also interesting that a risk model we proposed exhibited an excellent ability to predict the prognosis of TC. We expected to identify several clinically applicable diagnostic or prognostic biomarkers, and we provided a novel perspective of autophagy in TC.
Due to the indispensable function of autophagy, some researchers have evaluated several effects to reveal the clinical significance and corresponding molecular mechanisms in TC. However, its role in the tumorigenesis of TC remains controversial and unclear. Several scientific researchers have studied autophagy-related agents combined with other clinical treatment techniques in TC (31,32). These strategies regarding autophagy inhibition have achieved certain effects; however, there are still issues that have not been resolved. Hence, it is imperative to analyze the available data to indicate that the role of autophagy in the steps of TC onset, progression and prognosis. In the present study, a differential expression analysis was conducted among a series of ARGs to select the key mediators in the initiation of TC. In addition, we performed a functional analysis to reveal the unequivocal role of autophagy. These results revealed that in addition to autophagy, these differentially expressed genes were also involved in several proliferation-related pathways. Notably, the proliferative signaling pathways were the core of the interference. This conformed to the coordinating function of autophagy in the balance between cell growth and cell death (33–35). Not surprisingly, the DO analysis revealed that these differentially expressed ARGs also participated actively in various types of cancer. Further, the PPI network analysis indicated the existance of s interactions between these genes. BCL2, the first identified mammalian apoptotic regulator (36,37), is ranked as the core in the network. Moreover, BCL2-selective inhibitors have been approved for preclinical studies in anti-cancer treatment (38–40). The results of the functional enrichment analyses provided us with good reason to believe that these genes are common oncogenes, and that they operate as a whole to participate in the progression of the TC.
As the main purpose of this study was to reveal the clinical significance of autophagy in TC, we comprehensively analyzed the expression profiles, clinical parameters and prognostic values of these genes. As regards the 31 differentially expressed ARGs, 25 of these were significantly related to lymph node metastasis. Lymph node metastasis occurs frequently in patients with TC, and it directly affects the choice of treatment. Some previous studies have pointed to lymph node metastasis as an important indicator of TC invasion, and it was closely related with recurrence and distant metastasis (41–43), whether or not the lymph node metastasis had a great influence on the prognosis. Although there have been few studies reporting that several ARGs correlated with the lymph nodes, the majority of these research studies focused only on a signal gene (22,44,45). Our discovery that autophagy was associated with lymph node metastasis provided clues for the clinical management. However, further experimental validation, either in vivo or in vitro, is required in the future. These genes also displayed a significant correlation with the advanced tumor stage. These findings evidently show that these genes actively participate in the proliferation and invasion of TC. Therefore, these genes may be particularly useful for the assessment of the treatment response and progression surveillance.
We examined the associations between the ARG expression levels and survival using a univariate Cox regression model, and identified 6 genes as candidate prognostic biomarkers. We then submitted these genes to a multivariable Cox regression and developed a signature for predicting survival. The signature made it possible to separate the patients with TC into two groups with significantly different survivals. Three genes (ATG9B, BID and DNAJB1) included in the prognostic signature could act as clinically applicable indicators. Autophagy process is under strict control by a series of ARGs (46). Previously, it was demonstrated that the dysregulated autophagy level was closely related to tumor growth, survival and proliferation (47). Moreover, the deletion of several essential ARGs compromises the survival of tumor cells in vitro and in vivo (48,49). Therefore, the stable expression of ARGs is indispensable for inhibiting the onset of tumors. In this study, we found that several ARGs were dysregulated in the onset and progression of PTC. Therefore, these genes may play an important role in leading uncontrolled autophagy. Previously, Wang et al reported that the lack of ATG9B may facilitate the docking of both LC3 and p62 to initiate autophagy-associated degradation (50). Furthermore, this process may promote the apoptosis of hepatocytes in the initiation of hepatocellular carcinoma. Behrends et al first reported that DNAJB1 was involved in autophagy, which served as the bridge between WIPI2 and ATG2A. WIPI2 and ATG2A act with each other through DNAJB1. And depletion of WIPI2 led to reducing numbers of autophagosomes (51). In MCF-7 breast cancer cells, BID knockdown has been proved effective in the suppression of apoptosis and a shift of cell death towards autophagy (52). These findings displayed different roles of core ARGs in autophagy. Considering their expression profiles and clinical significance in PTC, these genes may also contribute to the progression of PTC via autophagy. ATG9B, BID and DNAJB1 all play multiple functions in autophagy; however, they have not been reported for their important clinical significance and molecular function in TC. These genes of less concern could provide novel insight into the clinical management and molecular mechanisms responsible for the development and progression of TC. Therefore, in this study, we used several bioinformatics methods to dig up the potential regulatory mechanism. These three genes are all significantly related to the methylation level of the CpG locus. Of note, the patients with alterations in these three genes had worse survivals when compared with the patients without alterations. Thus, we proposed the hypothesis that epigenetics and genomic alterations change the expression profiles of ARGs, thus influencing the prognosis of TC. In the present study, we also explored the associations between alterative splicing events and core ARGs. We also found the difference of alterative splicing events in each tumor samples were not huge. These findings hinted that alterative splicing events may not be the cause of significant clinical values of the three genes.
In this study, some limitations should be acknowledged. First, the underlying molecular mechanisms of the key ARGs in the TC pathogenesis are lacking, and additional experimental investigations should be conducted to uncover these. Second, although the prognostic signature we proposed was based on the expression profiles and clinical information obtained from the generally recognized TCGA database, it should be validated in other independent databases. Despite these limitations, this study highlighted the great clinical significance and potential molecular mechanisms of autophagy. These findings suggest that autophagy-targeted treatment may have a unique effect on TC, particularly in patients with lymph node metastases, and they provide insight into the complex biological functions of autophagy in TC.