This article is mentioned in:

Introduction

As the most common thyroid malignancy worldwide, papillary thyroid carcinoma (PTC) accounts for ~80% of all thyroid cancers (1). The incidence rate of PTC exhibits the most rapid increase of all cancers among women, and the second most rapid increase among men (2). Despite the advances in surgery, ~5% of patients with PTC experience recurrence within 5 years of the initial treatment (3). Understanding the molecular mechanism involved in the proliferation, apoptosis and invasion of PTC is extremely important for the development of more effective therapeutic strategies.

Recently, certain differentially-expressed genes (DEGs) have been reported to exhibit important roles in PTC, and identification of these may be useful in the investigation of the molecular mechanisms of PTC (4,5). Previous studies showed that ret proto-oncogene, neurotrophic tyrosine kinase receptor type 1 and v-raf murine sarcoma viral oncogene homolog B may be useful therapeutic targets for PTC (6–8). Ye et al found that Krüppel-like factor 17 may serve as a candidate tumor suppressor and a therapeutic target in PTC (9). Programmed cell death 4 was reported to exhibit an inhibitory role in the cell proliferation, malignant progression and invasion of PTC (10). microRNA (miRNA/miR)-199b-5p, miR-30a-3p and miR-146b-5p may be associated with PTC invasiveness (11). Although serious attempts have been made to find novel targets for gastric cancer treatment, at present, this knowledge is insufficient.

In the present study, DEGs between PTC patients and normal individuals were identified. Modules were then screened from the protein-protein interaction (PPI) network and the significant target genes were selected from the miRNA regulatory network. Through the identification of key genes, the possible molecular mechanism and potential therapeutic targets for PTC were investigated.

Materials and methods

Affymetrix microarray data

The gene expression profile, GSE53157, which was deposited by Pita et al (12), was obtained from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/). Gene expression profiling was based on the platform of GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array). The array consists of 54,675 probesets that can be used to detect the transcription level of 18,750 human genes. Only 10 chips, including 3 specimens of normal thyroid tissues and 7 specimens of well-differentiated thyroid carcinomas, were analyzed in the present study.

Identification of DEGs

The raw data were first preprocessed using the Affy package (13) in R language. Next, DEGs between normal thyroid tissues and well-differentiated thyroid carcinomas were analyzed by limma package in R (14). Fold-change (FC) of the expression of individual genes was also calculated for the differential expression test. DEGs with an adjusted P-value (Adj.P.Val) of <0.05 and |log FC| ≥1 were considered to be significant. Adj.P.Val was the result of multiple testing corrections using the Benjamini-Hochberg (HB) method (15).

Gene Ontology (GO) and pathway enrichment analysis of DEGs

The GO analysis has become a commonly used approach for functional studies of large-scale transcriptomic or genomic data (16). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (17) contains information of how molecules or genes are networked. The Database for Annotation Visualization and Integrated Discovery (DAVID) (18) was used to systematically extract biological meaning from large gene or protein lists. The GO function and KEGG pathway of DEGs were analyzed using DAVID 6.7 with FDR<0.05.

Construction of PPI network and screening of module

Genes associated with PTC were downloaded from the Malacards database (19). The downloaded genes and the identifided DEGs were then combined; the pooled dataset is referred to as PTC-associated genes in the present study. The Search Tool for the Retrieval of Interacting Genes (STRING) (20) database was used to retrieve the predicted interactions for the PTC-associated genes; version 9.1 of STRING covers 1,133 completely sequenced species. All associations obtained in STRING are provided with a confidence score, and each score represents a rough estimate of how likely a given association describes a functional linkage between two proteins (21). The PTC-associated genes with a confidence score of >0.4 were selected to construct the PPI network, using the Cytoscape software (22).

Visualizing complex networks and integrating these networks to any type of attribute data was allowed by Cytoscape (http://cytoscape.org/). The clusterMaker 1.11 (23) plugin in Cytoscape and the Markov cluster algorithm (24) were used to divide the PPI into modules (granularity parameter, 2). GO functional analysis of the genes in the modules was performed using the BinGo 2.44 plugin in Cytoscape (25) with a threshold of P<0.05 using the hypergeometric test.

Enrichment analysis of microRNA targets

The PTC-associated microRNAs were screened from genes that were obtained from the Malacards database (19). The identified DEGs were then submitted into a web-based gene set analysis toolkit (WebGestalt) (26), and miRNA-target gene enrichment analysis was performed. The enrichment significance of the predicted target genes in the DEGs was tested using hypergeometric distribution. The target genes with an adj.P.Val of <0.01 were considered to be significant, which was the result of multiple testing correction using the HB method (15). Next, miRNA regulatory networks were built using Cytoscape (25). Furthermore, the integrated network was constructed using the miRNA regulatory network and the PPI modules.

Results

GO and pathway enrichment analysis of DEGs

In total, 668 DEGs, including 262 upregulated genes and 406 downregulated genes, were selected. Results of GO analysis showed that the upregulated DEGs were significantly enriched in biological processes such as the regulation of protein kinase activity, the regulation of transferase activity and the induction of programmed cell death (Table I). Results of pathway analysis showed that the upregulated DEGs were significantly enriched in the p53 signaling pathway (Table I). However, there were no significant GO biological processes and pathways in the downregulated DEGs.

Table I. GO and pathways analysis of the upregulated DEGs.

Table I.

GO and pathways analysis of the upregulated DEGs.

Category	Term	Count	FDR
GOTERM_BP_FAT	GO:0045859-regulation of protein kinase activity	18	1.17×10−2
GOTERM_BP_FAT	GO:0045860-positive regulation of protein kinase activity	14	1.22×10−2
GOTERM_BP_FAT	GO:0043549-regulation of kinase activity	18	1.36×10−2
GOTERM_BP_FAT	GO:0033674-positive regulation of kinase activity	14	1.47×10−2
GOTERM_BP_FAT	GO:0051347-positive regulation of transferase activity	15	1.51×10−2
GOTERM_BP_FAT	GO:0012502-induction of programmed cell death	16	1.96×10−2
GOTERM_BP_FAT	GO:0010942-positive regulation of cell death	19	2.16×10−2
GOTERM_BP_FAT	GO:0006917-induction of apoptosis	16	2.35×10−2
GOTERM_BP_FAT	GO:0051338-regulation of transferase activity	19	2.54×10−2
GOTERM_BP_FAT	GO:0043068-positive regulation of programmed cell death	19	2.61×10−2
GOTERM_BP_FAT	GO:0043627-response to estrogen stimulus	9	2.93×10−2
GOTERM_BP_FAT	GO:0043065-positive regulation of apoptosis	18	4.20×10−2
KEGG_PATHWAY	hsa04115:p53 signaling pathway	9	3.46×10−2

[i] GO, Gene Ontology; DEG, differentially-expressed gene; BP, biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes; count, number of DEGs in each GO term; FDR, false discovery rate.

Module screening from the PPI network

A total of 355 genes asscoiated with PTC were downloaded from the Malacards database. These genes and the identified DEGs were combined, and 983 PTC-associated genes were obtained. The PPI network was constructed based on the predicted interactions of the 983 PTC-associated genes. A total of 1,157 unique PPI pairs and 413 nodes were included. The PPI network was divided into three modules (Fig. 1A–C). Tumor protein p53 (TP53) was the core of module 1, which included 55 unique PPI pairs and 31 nodes. In module 1, the upregulated genes of cyclin-dependent kinase inhibitor 1A and S100 calcium binding protein A6 (S100A6) were included in the Malacards database (Fig. 1A). V-raf-1 murine leukemia viral oncogene homolog 1 was the core of module 2, and a total of 45 unique PPI pairs and 26 nodes were included in module 2. In module 2, the upregulated genes of dual specificity phosphatase 5 and keratin 19 were included in the Malacards database (Fig. 1B). Epidermal growth factor receptor was the core of module 3, and a total of 39 unique PPI pairs and 21 nodes were included in this module. In module 3, the upregulated genes of met proto-oncogene (MET) and lectin galactoside-binding soluble 3 were included in the Malacards database (Fig. 1B). Modules 1, 2 and 3 were found to be significantly enriched (P<0.05) for one GO term each (Table II). The significant GO term in module 1 was the negative regulation of cellular metabolic process (P=8.16×10−6). The term response to chemical stimulus was the most significantly enriched function in module 2 (P=7.53×10−7). Notably, the transmembrane receptor protein tyrosine kinase signaling was the most significantly enriched pathway (P=6.50×10−11) in module 3.

Figure 1. Selected modules from the protein-protein interaction network: (A) Module 1, (B) module 2 and (C) module 3. Orange and blue nodes represent products of up- and downregulated DEGs, respectively. Grey nodes represent PTC-associated genes that were not included in the DEGs. Nodes with a purple edge represent PTC-associated genes that were collected in the Malacard database. Nodes with a green edge represent PTC-associated genes that were not collected in the Malacard database. DEG, differentially-expressed gene; PTC, papillary thyroid carcinoma.

Table II. GO terms of genes in three modules.

Table II.

GO terms of genes in three modules.

Module	GO term	Corr. P	Genes in test set
1	Negative regulation of cellular metabolic process	8.16×10−6	E2F1, CDKN1A, CDKN2A, PHB, TRIM27, TP53, MDM2, TP63, MDM4, BHLHE40, RUNX2, STAT3
2	Response to chemical stimulus	7.53×10−7	HNF1A, BRAF, PMAIP1, BCL2L1, SOD1, MAPK1, DUSP4, KRT19, KRAS, BAX, BCL2, MAPK3, SCNN1G, SCNN1B, DUSP6
3	Transmembrane receptor protein tyrosine kinase signaling pathway	6.50×10−11	PTPRJ, EGFR, RET, ERBB3, ERBB2, GAB1, PDGFRB, EGF, NRG1, GRB7

[i] GO, Gene Ontology; Corr. P, corrected P-value.

Enrichment analysis of microRNA target genes

The PTC-associated microRNAs acquired from the 355 genes were miR203, miR34b, miR221, miR222, miR146a, miR146b, miR21, miR181a1 and miR130b. Through miRNA-target gene enrichment analysis, a miRNA regulatory network consisting of 54 miRNAs and 210 DEGs was constructed (Fig. 2A). miR-221, miR-222, miR-181a, miR-203 and miR-130b, which were included in the 355 genes, were involved in the miRNA regulatory network. In addition, cyclin-dependent kinase inhibitor 1C (CDKN1C) was the target of miRNA-221 and miRNA-222; v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog, metallophosphoesterase domain containing 2 and peroxisome proliferator-activated receptor γ were the targets of miRNA-130b; aryl hydrocarbon receptor was the target of miRNA-203; and dual specificity phosphatase 5 was the target of miRNA-203 and miRNA-181a (Fig. 2B). By combining the miRNA regulatory network with the PPI network modules, the integrated network was constructed (Fig. 3A–C). The upregulated gene basic helix-loop-helix family member e40 in module 1 was the target of miR-181a (Fig. 3A); and the downregulated gene reticulon 1 in module 2 was the target of miR-130b (Fig. 3B).

Figure 2. miRNA regulatory network. (A) miRNA regulatory network constructed with DEGs and their regulatory miRNAs. (B) miRNA regulatory network constructed with DEGs and their regulatory miRNAs. Diamond-shaped nodes represent miRNA and circular nodes represent DEGs. Orange and blue nodes represent products of up- and downregulated DEGs, respectively. Grey nodes represent PTC-associated genes that were not included in the DEGs. Nodes with a purple edge represent PTC-associated genes thatwere collected in the Malacard database; and nodes with a green edge represent PTC-associated genes that were not collected in the Malacard database. miRNA, microRNA; DEG, differentially-expressed gene; PTC, papillary thyroid carcinoma.

Figure 3. Integrated regulatory network. (A) Integrated network of module 1 and miRNA regulatory network. (B) Integrated network of module 2 and miRNA regulatory network. (C) Integrated network of module 3 and miRNA regulatory network. Diamond-shaped nodes represent miRNA and circular nodes represent DEGs. Orange and blue nodes represent products of up- and downregulated DEGs, respectively. Grey nodes represent PTC-associated genes that were not DEGs. Nodes with a purple edge represent PTC-associated genes that were collected in the Malacard database; and nodes with a green edge represent PTC-associated genes thatwere not collected in the Malacard database. miRNA, microRNA; DEG, differentially-expressed gene; PTC, papillary thyroid carcinoma.

Discussion

PTC accounts for nearly 80% of human thyroid cancers, worldwide (27). Now, the burden of thyroid disease in the global population is enormous (28). Thus, the potential use of therapeutic targets appears to be the most promising area of research. In the present study, a bioinformatics approach was used to predict the potential therapeutic targets for PTC. The study identified 668 DEGs between normal thyroid tissues and well-differentiated thyroid carcinomas, among which, 406 genes were downregulated and 262 were upregulated. By constructing the PPI and miRNA regulatory networks, key genes were found, including S100A6, MET and CDKN1C.

S100A6 encodes a member of the S100A family of calcium binding proteins, and it is reported to be involved in the regulation of a wide range of cellular processes, including cell proliferation, differentiation and apoptosis (29,30). Furthermore, S100A6 has been shown to be upregulated in numerous cancers, including thyroid carcinoma (31,32). Li et al showed that increased S100A6 expression promoted the proliferation and migration of cells in human hepatocellular carcinoma (HCC), and that S100A6 may function as a potential therapeutic target in HCC (33). Komatsu et al showed that S100A6 was involved in the invasion and metastasis processes of human colorectal adenocarcinomas (34). S100A6 has been reported to be overexpressed in PTC and may contribute to certain events in this disease (35,36). In the present study, S100A6 was overexpressed in PTC, and the PPI network showed that it had direct interaction with TP53, which was associated with PTC (37,38). Overall, S100A6 may be a therapeutic target in PTC.

MET encodes the hepatocyte growth factor receptor and could regulate a number of physiological processes, including proliferation, scattering (cell dissociation and motility), morphogenesis and survival (39). MET could trigger tumor growth, angiogenesis and metastasis (40). Rong et al found that MET was overexpressed in a number of tumors, and suggested that MET activation contributed to tumor progression (41). Fujita and Sugano found that MET was overexpressed in primary colorectal cancer and liver metastases, and suggested that MET played an important role in the development of colorectal cancer liver metastases (42). In one study, Di Renzo et al showed that the overexpression of MET in human thyroid carcinomas may contribute to the progression of thyroid tumors (43), while in a second study, it was suggested that the upregulation of MET may add a selective growth advantage to differentiated ovarian cancers (44). In the present study, MET was upregulated in PTC, and GO analysis showed that it was involved in the biological process of transferase activity. Combined with the aforementioned studies, this data shows that MET may be a therapeutic target in PTC.

The putative tumor suppressor CDKN1C is a negative regulator of cell proliferation and mutations. Soejima et al showed that CDKN1C was downregulated in esophageal cancer and contributed to the tumor (45), and Hoffmann et al found that numerous advanced urothelial cancers displayed the downregulation of CDKN1C (46). Larson et al found that CDKN1C was a candidate tumor suppressor in breast cancer (47), while Algar et al found that it was a tumor suppressor in rhabdoid tumors (48). However, there have been no studies with regard to the association between CDKN1C and PTC. In the present study, CDKN1C was downregulated in PTC, and the miRNA regulatory network showed that CDKN1C was the target of miRNA221 and miRNA222. One previous study showed that miRNA221 and miRNA222 were upregulated in PTC, and that they contributed to the procession of the disease (49). Furthermore, Fornari et al indicated that miR-221 controlled the expression of CDKN1C in human hepatocellular carcinoma (50). Thus, CDKN1C may be a potential therapeutic target in PTC.

In conclusion, S100A6, MET and CDKN1C may exhibit key roles in the progression and development of PTC, and may be used as specific therapeutic targets in the treatment of PTC. However, further experiments are required to confirm these results.

Acknowledgements

The authors wish to express gratitude to Fenghe (Shanghai) Information Technology Co., Ltd. (Shanhai, China) for ideas and assistance that provided a valuable added dimension to the present study.

References