This article is mentioned in:

Introduction

Thyroid cancer is the most common type of malignancy in the head and neck. It can be divided into four pathological types: Papillary thyroid carcinoma (PTC), follicular carcinoma, medullary carcinoma and undifferentiated carcinoma. PTC is the most common of these and accounts for ~80% of cases (1). The incidence of PTC in the USA has increased in recent years (2). Nevertheless, the prognosis for PTC is good and overall survival is ~90% (3).

Alterations in the RET/PTC-RAS-BRAF signaling pathway are characteristic of PTC (4,5). Activating BRAF and ret/PTC1 mutations may be an important step in the development of this disease (6–10). The etiology of PTC remains poorly understood. Therefore the identification of diagnostic and prognostic biomarkers is required. The detection of thyroid nodules is facilitated by high-resolution sonography. However, the detection of the nodules as either benign or malignant using fine-needle aspiration biopsy, real-time sonographic elastography, ultrasound elastography or Micro-Pure imaging, is problematic (11).

Microarray technology may be an effective tool for investigating the underlying regulatory network involved in PTC (12–15). In addition, He et al (16) have elucidated the role of microRNAs in the predisposition to, and development of PTC. Despite this, a greater understanding of PTC is required in order to improve clinical diagnosis and treatment, for example by identifying early biomarkers.

In the present study, gene expression profiles of thyroid tissue from patients with PTC were compared with those from healthy controls, in order to identify differentially expressed genes (DEGs). The regulatory associations for the DEGs were retrieved and potential key genes were identified. Pathway enrichment and biological process analyses were conducted to investigate the roles of these genes in the development of PTC.

Materials and methods

Microarray data

Microarray data set, GSE3678, consisting of 14 thyroid tissue samples, was downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). This included seven thyroid tissue samples from patients with PTC and seven from healthy people. Raw data were obtained from Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Inc., Santa Clara, CA, USA).

Screening of DEGs

Raw data were processed using software package GEOquery in R and differential expression analysis was conducted using limma (17) in R. The Beyer-Hardwick (BH) method was adopted for multiple corrections (18) in order to identify DEGs between the disease and the control groups. A fold change (FC) of >2 (| logFC |> 2) and P<0.05 were set as the cutoff values.

Establishment of a gene regulatory network

Gene regulation analysis was conducted using TRANSFAC® and TRED (19). A total of 5,558 regulatory associations were identified in TRANSFAC and TRED and extracted to a single database. The 171 DEGs were entered into this database in order to set up a regulatory network, and the corresponding schematic diagram of 104 DEGs was generated using Cytoscape version 2.8 (20).

Pathway enrichment analysis, biological process analysis and gene ontology (GO) annotation

Pathway enrichment analysis and biological process analysis with GO annotation were conducted for genes included in the regulatory network, using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (21) in order to determine their involvement in the development of PTC.

Results

DEGs

Differential expression analysis identified 171 genes as DEGs.

Gene regulatory network

The regulatory associations for the 171 DEGs were identified from the pre-established database. In total, 111 associations between transcription factors (TFs) and target genes among 104 DEGs were determined by the comparison between 5,558 regulatory associations and 171 DEGs, using TRANSFAC and TRED (Fig. 1).

Figure 1 Regulatory network for DEGs. Squares indicate transcription factors while circles represent target genes. DEGs are shown in yellow and non-DEGs are shown in pink. DEGs; differentially expressed genes.

Among the 104 DEGs, three TFs were identified: Trefoil factor 3 (TFF3), cut-like homeobox 2 (CUX2) and forkhead box protein A2 (FOXA2). The target genes of these DEGs are listed in Table I.

Table I Differentially expressed transcription factors and their target genes.

Table I

Differentially expressed transcription factors and their target genes.

TF gene	Target gene
TFF3	IFNB1
CUX2	CYBB
FOXA2	APOB
FOXA2	F2
FOXA2	IGFBP1
FOXA2	PROC
FOXA2	TTF1
FOXA2	MSH2

[i] TF, transcription factor; CUX2, cut-like homeobox 2; FOXA2, forkhead box protein A2; IFNB1, interferon β 1; CYBB, cytochrome b-245, β polypeptide; APOB, apolipoprotein B; IGFBP1, insulin-like growth factor binding protein 1; TTF1, thyroid transcription factor 1; MSH2, mutS homolog 2.

Pathway enrichment analysis, biological process analysis and GO annotation

Pathway enrichment analysis was conducted for the 104 DEGs. It identified six pathways (Table II); the most significant of these was in cancer. These results suggest that certain DEGs participate in pathways involved in PTC.

Table II Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis results for the 104 DEGs.

Table II

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis results for the 104 DEGs.

Term ID	Pathway name	Count	P-value	FDR
hsa05200	Pathways in cancer	26	4.72E-15	5.04E-12
hsa05221	Acute myeloid leukemia	10	1.62E-08	1.71E-05
hsa05220	Chronic myeloid leukemia	9	2.12E-06	0.002235
hsa05222	Small cell lung cancer	9	5.04E-06	0.005318
hsa05212	Pancreatic cancer	8	1.80E-05	0.019025
hsa05216	Thyroid cancer	5	3.34E-04	0.351572

[i] DEGs, differentially expressed genes; FDR, false discovery rate.

Biological process analysis was also conducted for the 104 DEGs. Certain DEGs were associated with the positive regulation of gene expression, gene transcription and metabolic processes (Table III).

Table III GO biological process analysis results for the 104 DEGs.

Table III

GO biological process analysis results for the 104 DEGs.

GO Term	Description	Count	P-value	FDR
GO:0010628	Positive regulation of gene expression	43	1.43E-31	2.40E-28
GO:0045944	Positive regulation of transcription from RNA polymerase II promoter	37	2.78E-31	4.67E-28
GO:0045893	Positive regulation of transcription, DNA-dependent	40	3.57E-31	6.00E-28
GO:0051254	Positive regulation of RNA metabolic process	40	4.92E-31	8.27E-28
GO:0045941	Positive regulation of transcription	42	7.54E-31	1.27E-27
GO:0010557	Positive regulation of macromolecule biosynthetic process	43	1.77E-29	2.98E-26
GO:0045935	Positive regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolic process	42	4.21E-29	7.07E-26
GO:0051173	Positive regulation of nitrogen compound metabolic process	42	1.46E-28	2.46E-25
GO:0009891	Positive regulation of biosynthetic process	43	2.06E-28	3.46E-25
GO:0006357	Regulation of transcription from RNA polymerase II promoter	43	1.25E-27	2.10E-24

[i] DEGs, differentially expressed genes; GO, gene ontology; FDR, false discovery rate.

Discussion

In the present study, microarray data for thyroid tissue samples from PTC patients and from healthy controls were compared. A total of 171 DEGs were identified in these samples. Gene regulation analysis revealed 111 associations among 104 of these DEGs. Three differentially expressed TFs were identified: TFF3, CUX2 and FOXA2. Pathway enrichment and biological process analyses were conducted for the 104 DEGs in order to investigate their involvement in the development of PTC. Pathway enrichment analysis indicated that these genes are associated with colorectal, pancreatic, thyroid and lung cancers, and leukemia (Fig. 2). Biological process analysis suggested that these genes exhibit corresponding functions through positive regulation of gene expression, transcription and metabolic processes.

Figure 2 Significant pathways for genes in the regulatory network: Pathways in cancer.

Underexpression of TFF3 in PTC has been confirmed by previous studies (22,23) and it is proposed as a biomarker at the RNA level (24,25). Its target gene is interferon β 1 (Table I), which exhibits anticancer activity (26). CUX2 (Table I) contains three CUT domains, which are three internal CUT repeats, and a homeodomain, which are DNA-binding motifs. CUX2 is important in the negative regulation of transcription and has been linked to tumor progression (27). FOXA2 is a member of the forkhead box-O family (Table I), which is involved in the regulation of a series of cellular events, such as differentiation, DNA repair, cell cycle arrest and apoptosis (28,29). It is a potential oncogene in anaplastic thyroid carcinoma (30). Kim et al (31) confirmed the decreased expression of FOXA2 in PTC cells using reverse transcription-quantitative polymerase chain reaction. Akagi et al (32) also found a reduced expression of FOXA2 in PTC cells compared with that in healthy thyroid cells. The enforced expression of a range of genes, including FOXA2, inhibits PTC cell growth (32). Furthermore, Akagi et al (32) reported that the CpG island in the promoter region of FOXA2 was aberrantly methylated and that treatment with 5-aza-2-deoxycytidine induced the expression of FOXA2 in PTC cells. Insulin-like growth factor binding protein 1 (IGFBP1) is a target gene of FOXA2 (Table I). Yashiro et al (33) demonstrated that IGFBP1 activity is significantly higher in cells from patients with PTC, therefore it may be involved in PTC cell growth regulation.

Pathway enrichment analysis indicated that the top pathway that the DEGs were found to be associated with was pathways in cancer, involving 26 DEGs. Fibronectin 1 is a glycoprotein present in a dimeric or multimeric form, which is present at the cell surface and in the extracellular matrix. It is involved in cell adhesion and migration, and is therefore a potential drug target with which to treat cancer, by blocking the processes involved in metastasis (34). Its overexpression in PTC at the mRNA level has been demonstrated by Takano et al (35). Cyclin-dependent kinase inhibitor 2B forms a complex with cyclin-dependent kinase 4 or cyclin-dependent kinase 6, and functions as a cell growth regulator (36). Arnaldi et al (37) revealed that it is underexpressed in malignant thyroid tissue. Peroxisome proliferator-activated receptor γ 1 (PPARG1) fuses with paired box 8 and has been confirmed as an oncogene in human thyroid carcinoma (38). In the present study, peroxisome proliferator-activated receptor γ coactivator 1-α, a target gene of PPARG1, was found to be differentially expressed, which was in accordance with the findings of Antico Arciuch et al (39).

Biological process analysis suggested that the DEGs identified in the present study were linked to gene expression, gene transcription, and the positive regulation of RNA and nitrogen compound metabolic processes (Table III). Therefore, a variety of regulatory mechanisms were involved in the development of PTC. Cbp/p300-interacting transactivator, with Glu/Asp-rich carboxy-terminal domain 1 (CITED1) and CITED2 were transcriptional coactivators and were identified as DEGs. Prasad et al (14) indicated that CITED1 was upregulated in PTC cells compared with cells from healthy thyroids. Therefore it may serve as a diagnostic marker for PTC.

In conclusion, 171 DEGs in PTC were identified and the molecular mechanisms of 104 DEGs were investigated following the creation of a regulatory network. Three TFs associated with PTC were identified. Pathway enrichment and biological process analyses of the genes in the network confirmed their involvement in the development of PTC. The results of the present study were beneficial for future research into the genes and regulatory mechanisms associated with the pathogenesis of PTC.

Acknowledgements

This study was supported by a grant from Capital Clinical Feature Application Research (Z141107002514102).

References