This article is mentioned in:

Introduction

Thyroid cancer is the most common type of cancer of the endocrine system, with cases increasing worldwide (1,2). Thyroid cancer may be classified into numerous types according to the histopathological characteristics. Medullary thyroid carcinoma (MTC) is a form of thyroid cancer which originates from the parafollicular cells of the thyroid (3). It is the third most common type of thyroid cancer and accounts for ~3% of all thyroid cancer cases. Approximately 1 in 4 of MTC cases are caused by mutations in the rearranged during transfection (RET) proto-oncogene (4). The majority of MTC cases are sporadic, presenting with metastatic disease at diagnosis (5). Nearly all patients with distant metastases succumb to this disease (6).

Presently, RET mutation have been suggested to be an indicator of the poor prognosis of MTC; however, this is not sufficient for understanding the underlying molecular mechanisms of MTC tumourigenesis. Soh et al (7) demonstrated that vascular endothelial growth factor receptor 2 was involved in the pathogenesis of MTC via promotion of pro-invasive and pro-angiogenic phenotypes. Additionally, dysregulation of the Dickkopf/Wnt signaling pathway inhibitor 4 has been identified in MTC (8). A previous study revealed that aberrant expression levels of microRNAs (miRNAs) have a potential role in tumourigenesis (9), which may provide novel insight in MTC research. Notably, increasing evidence has supported the important role of miRNAs in cancers including thyroid cancer (10,11). For example, He et al (12) reported that three miRNAs, including miR-221, −222, and −146, are overexpressed in papillary thyroid cancer. Furthermore, miR-197 and −346 are significantly overexpressed in follicular thyroid cancers (13). Although great advances have been made in understanding the functions of miRNAs in thyroid cancers, the underlying molecular mechanisms of this disease remain to be elucidated.

The present study aimed to use GSE40807 miRNA microarray data provided by Lassalle et al (14) to identify differentially expressed miRNAs (DEMs) between human MTC and healthy control tissues. Subsequently, transcription factor (TF)-miRNA and miRNA-target gene regulatory networks were constructed. Finally, the target genes of DEMs were performed functional enrichment analyses to predict their potential functions that may be associated with MTC. To the best of our knowledge, this is the first time that the dataset of GSE40807 was analyzed.

Materials and methods

Data source

The GSE40807 miRNA microarray dataset was downloaded from the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/) database in the National Center for Biotechnology Information based on the Agilent-019118 Human miRNA Microarray 2.0 G4470B platform (Agilent Technologies, Inc., Santa Clara, CA, USA). The dataset included 14 pairs of miRNA microarrays from human MTC and adjacent healthy tissues.

Data preprocessing and DEM analysis

The original data were firstly converted into identifiable expression form using a Linear Models for Microarray Data (limma) package in R language (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (15). Following this, background correction and quartile data normalization were performed using a robust multiarray average algorithm affy package in R (http://www.bioconductor.org/packages/release/bioc/html/affy.html) (16).

TF-miRNA regulatory association pair prediction

The TF-miRNA regulatory database (TransmiR; cmbi.bjmu.edu.cn/transmir) (17) is a valuable resource for the study of TF-miRNA regulation, which provides an interface for easy retrieval of TF-miRNA regulatory pairs by searching for a miRNA or a TF. Currently, TransmiR has curated 735 entries, which includes 201 TFs, 209 miRNAs and 16 organisms from 268 publications. The present study inputted the obtained DEMs into the database and extracted regulatory association pairs between DEMs and TFs.

miRNA-target gene regulatory association pair prediction and TF-miRNA-target gene regulatory network construction

The starBase v2.0 (starbase.sysu.edu.cn/) database (18) provides certain miRNA-target regulatory association pairs which are verified by experiments and predicted by five algorithms including TargetScan (19), miRanda (20), Pictar2 (21), PITA (22) and RNA22 (23). In the present study, miRNA-target gene regulatory association pairs verified by ≥ one experiment and predicted by ≥ three algorithms were selected for construction of the regulatory network.

Based on the predicted TF-miRNA and miRNA-target gene regulatory association pairs, the TF-miRNA-target gene regulatory network was constructed using Cytoscape software version 3.2.0 (www.cytoscape.org/) (24). From the network, TFs, miRNAs and target genes that had higher connective degrees (hub nodes) were extracted. Hub nodes are small numbers of nodes with numerous interaction partners, which serve important roles in the network (25). Thus, these TFs, miRNAs and target genes might serve roles in MTC.

Functional enrichment analyses

clusterProfiler software (Bioconductor version 3.1; bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (26) is a package used for gene classification and enrichment analysis. The Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.ad.jp/kegg/) (27) is a database of biological systems that collects genomic, chemical and systemic functional information. To analyzed the potential biological functions of DEMs, KEGG pathway enrichment analysis was performed for the target genes of the obtained DEMs based on the clusterProfiler package. P<0.01 was set as the threshold value.

Statistical analysis

DEMs between MTC and healthy tissues were identified using the limma (15) package (Bioconductor version 3.1). Student's t-test in the limma package was used to compare DEM values, and fold changes (FCs) were calculated. miRNAs with P<0.01 and |log2FC|≥1 were selected as DEMs. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of DEMs

A total of 32 DEMs were identified between MTC and healthy tissues. Among these DEMs, 15 were upregulated and 17 were downregulated (Table I). Hierarchical clustering analysis of these DEMs and samples are presented in Fig. 1.

Figure 1. Heatmap of differentially expressed miRs. The upper color bar represents sample classes; red represents disease samples and green represents healthy controls. miR, microRNA.

Table I. Up- and downregulated miRs.

Table I.

Up- and downregulated miRs.

miR	log2FC	P
hsa-miR-9-5p	1.771272789	9.96E-05
hsa-miR-149-5p	−1.285205715	0.000830471
hsa-miR-708-5p	−1.14714938	0.003865191
hsa-miR-335-5p	2.093096976	0.003949694
hsa-miR-592	1.802025779	0.004011795
hsa-miR-875-5p	−1.148677692	0.004350708
hsa-miR-455-5p	−1.191505453	0.004403134
hsa-miR-590-5p	1.642054145	0.004446886
hsa-miR-96-5p	2.177657366	0.005246821
hsa-miR-584-5p	−1.157995671	0.005083235
hsa-miR-922	1.022980071	0.005377546
hsa-miR-1	1.816907981	0.005578293
hsa-miR-296-5p	−1.335624531	0.00595956
hsa-miR-634	−1.222619212	0.006085225
hsa-miR-224-5p	1.794623266	0.006202874
hsa-miR-185-5p	1.573515208	0.006464766
hsa-miR-628-5p	1.035790796	0.007467156
hsa-miR-924	−1.075200672	0.007467315
hsa-miR-154-5p	1.230953928	0.007599707
hsa-miR-625-5p	1.081907512	0.009088003
hsa-miR-145-3p	1.656652903	0.00054755
hsa-miR-195-3p	−1.337668541	0.001551485
hsa-let-7f-1-3p	−1.15398965	0.001765947
hsa-miR-515-3p	−1.033833078	0.002354703
hsa-miR-9-3p	1.917564249	0.004292166
hsa-miR-7-2-3p	−1.106382457	0.005054543
hsa-miR-143-3p	1.498076967	0.005687211
hsa-miR-887-3p	−1.372273521	0.005784544
hsa-miR-512-3p	−1.149225488	0.006533263
hsa-miR-1237-3p	−1.009014763	0.006904609
hsa-miR-371a-3p	−1.099403781	0.007869777
hsa-let-7b-3p	−1.362383564	0.008104953

[i] log₂FC>0, upregulation; log₂FC<0, downregulation. miR; microRNA.

TF-miRNA regulatory relationship pairs

From TransmiR, 54 TF-miRNA regulatory association pairs were extracted, including 33 TFs and 10 miRNAs. Among these miRNAs, hsa-miR-9-5p was regulated by nine TFs, including nuclear factor of κ light polypeptide gene enhancer in B-cells 1 (NF-κB1), and interleukin 1β (IL-1β), and hsa-miR-1 was regulated by eight TFs, including CCAAT/enhancer binding protein α. Additionally, TFs of NF-κB1 regulated four DEMs, including hsa-miR-9-5p and −3p, and v-myc avian myelocytomatosis viral oncogene homolog (MYC) regulated three DEMs including hsa-miR-195-3p, hsa-let-7b-3p and hsa-let-7f-1-3p.

miRNA-target gene regulatory association pairs and TF-miRNA-target gene regulatory network construction

From starBase, 1654 miRNA-target gene regulatory association pairs were obtained, including 12 DEMs and 1338 target genes. Among the 12 DEMs, hsa-miR-1, hsa-miR-9-5p, hsa-miR-96-5p and hsa-miR-590-5p had the top four highest connective degrees (feature miRNAs). Additionally, seven target genes that were regulated by at least four DEMs were identified (Table II).

Table II. Feature miRNAs and genes in the miR-target gene regulatory network.

Table II.

Feature miRNAs and genes in the miR-target gene regulatory network.

Node	Number
hsa-miR-9-5p	326
hsa-miR-96-5p	297
hsa-miR-1	246
hsa-miR-590-5p	146
CRIM1	5
KIF1B	4
NR4A3	4
RNF111	4
TNPO1	4
FNDC3B	4
BCL11A	4

[i] miR, microRNA; CRIM1, cysteine rich transmembrane bone morphogenic protein regulator 1; KIF1B, kinesin family member 1B; NR4A3, nuclear receptor subfamily 4 group A member 3; RNF111, ring finger protein 111; TNPO1, transportin-1; FNDC3B, fibronectin type III domain containing 3B; BCL11A, B-cell lymphoma/leukemia 11A.

Furthermore, based on the constructed miRNA-target gene and TF-miRNA regulatory networks, a TF-miRNA- target gene regulatory network was constructed using Cytoscape software. In the network, there were 1654 miRNA-target gene and 54 TF-miRNA regulatory relationship pairs (Fig. 2).

Figure 2. TF-DEM-Target comprehensive regulatory network. Red squares represent microRNAs, orange triangles represent TF genes and blue circles represent non-TF genes. TF, transcription factor; DEM, differentially expressed microRNA.

Functional enrichment analyses

Among the 12 DEMs in the miRNA-target gene regulatory network, the target genes of hsa-miR-1, hsa-miR-9-5p, hsa-miR-96-5p and hsa-miR-590-5p were demonstrated to be enriched in the KEGG pathways, including the mitogen activated protein kinase (MAPK) signaling pathway, pathways in cancer, and during focal adhesion (Table III).

Table III. Enriched signaling pathways involving differentially expressed miRs.

Table III.

Enriched signaling pathways involving differentially expressed miRs.

miR	Pathway description
hsa-miR-1	Neurotrophin signaling pathway
	Renal cell carcinoma
	Axon guidance
	MAPK signaling pathway
	Pathways in cancer
hsa-miR-9-5p	Neurotrophin signaling pathway
	Bacterial invasion of epithelial cells
	Focal adhesion
	MAPK signaling pathway
	Endocytosis
hsa-miR-96-5p	GnRH signaling pathway
	ErbB signaling pathway
	Prostate cancer
	Neurotrophin signaling pathway
	Axon guidance
hsa-miR-590-5p	MAPK signaling pathway

[i] miR, microRNA; MAPK, mitogen activated protein kinase; GnRH, gonadotrophin-releasing hormone; ErbB, receptor tyrosine-protein kinase ERBB2.

Discussion

Patients with progressive MTC have limited treatment options (28). Thus, understanding the underlying molecular mechanism of carcinogenesis may facilitate diagnosis and therapy options of this disease. In the present study, 15 upregulated and 17 downregulated DEMs were identified. In the constructed TF-miRNA regulatory network, hsa-miR-9-5p was regulated by 9 TFs and hsa-miR-1 was regulated by 8 TFs. The TFs of NF-κB1 and MYC regulated 4 and 3 DEMs, respectively. Additionally, the above two miRNAs served key roles in the miRNA-target gene regulatory network. Their target genes were primarily enriched in the MAPK signaling pathway and during focal adhesion. These miRNAs and signaling pathways may be important biomarkers for MTC diagnosis and treatment.

In the miRNA-target gene regulatory network, hsa-miR-1 was upregulated, and its target genes, including MAPK1, were enriched in numerous signaling pathways associated with cancer, including MAPK. MAPKs are a family of protein kinases whose functions are conserved during evolution from unicellular organisms (29). The MAPK signaling pathway consists of numerous key signaling components and phosphorylation events which control multiple fundamental cell processes including proliferation, differentiation and apoptosis (30). This signaling pathway has been frequently identified in activated in human cancers, which leads to malignant phenotypes including autonomous cell proliferation (31). Notably, Zatelli et al (32) demonstrated that the growth of the TT MTC cell line depends on activation of the MAPK signaling pathway, which suggests its role in MTC. In addition, MAPK signaling is important in regulating cytokine signaling pathways (33). Cytokines are released in response to inflammation and immunity, and have important roles in cancer development and progression (34). It has been reported that undifferentiated thyroid cancer cells secrete cytokines (35). Taken together, the MAPK signaling pathway may serve important roles in MTC via hsa-miR-1 and its target gene MAPK1.

In addition to hsa-miR-1, hsa-miR-9-5p upregulated in the miRNA-target gene regulatory network, and was regulated by 9 TFs in the TF-miRNA regulatory network, including NF-κB1 and IL-1ß. NF-κB is a transcription regulator activated by various intra- and extracellular stimuli. A previous study demonstrated that NF-κB1 regulates the expression of genes involved in numerous processes, including proliferation and apoptosis (36). Inappropriate activation of NF-κB has been associated with numerous inflammatory diseases, whereas persistent inhibition of NF-κB may lead to delayed cell growth (37). NF-κB has been demonstrated to be associated with the development of colorectal (38), breast (39), bladder (40), prostate (41) and advanced thyroid (36) cancers. On the other hand, IL-1β, a member of the IL-1 cytokine family, is an important mediator of the inflammatory response. A previous study reported that inflammation is a critical component of tumor progression (42). Zeki et al (43) suggested that IL-1 regulates G1 cell cycle progression and arrest in papillary thyroid carcinoma cells. Therefore, NF-κB1, IL-1ß and their regulated DEM hsa-miR-9-5p may serve important roles in MTC progression.

Additionally, TFs of MYC were demonstrated to regulate 3 DEMs including hsa-miR-195-3p, hsa-let-7b-3p and hsa-let-7f-1-3p. MYC is a multifunctional, nuclear phosphoprotein which serves roles in cell cycle progression, apoptosis and cellular transformation (44). The MYC gene has been widely implicated in numerous human cancers (45,46). Khosla et al (47) revealed that MYC mRNA expression levels increased in apoptotic TT cells, suggesting its role in MTC. Its regulated miRNA hsa-miR-195-3p has been identified to be abnormally expressed in a variety of cancers. For example, levels were upregulated in breast cancer and downregulated in gastric, hepatocellular and bladder cancers (48–50). The roles of hsa-let-7b-3p and hsa-let-7f-1-3p in cancer remain to be elucidated; thus, it was hypothesized that these miRNAs may be involved in MTC via regulation from the MYC gene. Taken together, MYC and its regulated DEMs, hsa-miR-195-3p, hsa-let-7b-3p and hsa-let-7f-1-3p, may serve important roles in the development of MTC.

The present study identified numerous key miRNAs and TFs that may be associated with MTC using comprehensive bioinformatics methods. However, no experiments with tissues or cells were performed to validate the expression levels of these miRNAs and TFs; a key limitation of this study. Additionally, there were only 14 pairs of miRNA microarrays in the dataset. Further studies with experimental validations and more samples are required to validate these observations.

In conclusion, the results of the present study indicated that the DEMs hsa-miR-1, hsa-miR-9-5p and hsa-miR-195-3p may have the potential to be used as diagnostic and therapeutic targets of MTC. Additionally, hsa-miR-1 and its target gene MAPK1 may serve a role in MTC, involving in MAPK signaling pathway. Additionally, TFs of IL-1ß and MYC may be implicated in the development of MTC.

Acknowledgements

The present study was supported by the Beijing Municipal Science & Technology Commission (grant no. Z141107002514003), the Beijing Municipal Administration of Hospitals Clinical Medicine Development of Special Funding Support (grant no. XMLX201311) and the National Natural Science Foundation of China (grant no. 81473499).

References